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Abstract: Several problems arising in Economics and Finance are analyzed using concepts and quantitative 
methods from Physics. The disertation is organized as follows: 

In the first chapter, it is argued that in a closed economic system, money is conserved. Thus, by analogy 
with energy, the equilibrium probability distribution of money must follow the exponential Boltzmann-Gibbs 
law characterized by an effective temperature equal to the average amount of money per economic agent. The 
emergence of Boltzmann-Gibbs distribution is demonstrated through computer simulations of economic mod- 
els. A thermal machine which extracts a monetary profit can be constructed between two economic systems 
with different temperatures. The role of debt and models with broken time-reversal symmetry for which the 
Boltzmann-Gibbs law does not hold, are discussed. 

In the second chapter, using data from several sources, it is found that the distribution of income is described 
for the great majority of population by an exponential distribution, whereas the high-end tail follows a power 
law. From the individual income distribution, the probability distribution of income for families with two earners 
is derived and it is shown that it also agrees well with the data. Data on wealth is presentend and it is found 
that the distribution of wealth has a structure similar to the distribution of income. The Lorenz curve and Gini 
coefficient were calculated and are shown to be in good agreement with both income and wealth data sets. 

In the third chapter, the stock-market fluctuations at different time scales are investigated. A model where 
stock-price dynamics is governed by a geometrical (multiplicative) Brownian motion with stochastic variance 
is proposed. The corresponding Fokker-Planck equation can be solved exactly. Integrating out the variance, an 
analytic formula for the time-dependent probability distribution of stock price changes (returns) is found. The 
formula is in excellent agreement with the Dow-Jones index for the time lags from 1 to 250 trading days. For 
time lags longer than the relaxation time of variance, the probability distribution can be expressed in a scaling 
form using a Bessel function. The Dow-Jones data follow the scaling function for seven orders of magnitude. 



Contents 



I Acknowledgements! 



G^DijtriMti^nof wealthi 



H. Other distributions for income and wealthi 



I. Conclusions! 



Statistical mechanics of monevi 



^^Boitzman^^ibb^isMbutjo3 

^^^lenna^nachina 

E. Models with debt! 



^^oltzm^m^2uation|_ 



G^Non-Boltzm^n^^ 



2 
2 
2 
3 
4 
4 
5 
6 



En: 



H. Nonlinear Boltzmann equation vs. linear master equation! / 



! I. Conclusions! 
! II. Distribution of income and wealthi 



Distribution of stock-price fluctuation^ 



B. The model! 



CSoluti^noftheFok^ equation! 
D^^ath4ntegral_SolutT^ 

F. Asymptotic behavior for long time d 



G^_Asvmptotic be haviorforlargejo g-return c3 
H. Comp Miggfl-W ith Dow- Jones time series! 
I. Conclusions! 



15 
16 
17 

18 
18 
18 
19 
20 
20 
21 
23 
24 
25 



! IV. Path-integral solution of the Cox-IngersoU-Ross/Feller model! 26 



^^^^^^^^^^^^^^^^^^^^^^als] 

^^Exgonentia^jstributio^n 



8 ! Referenced 



29 



^^owe^a^ai^n^|Bose]^ondensatioi^ 



^^Geogra^hkaWariation^i^nco^ 
^^Distributioi^nncom^fo^fai^^ 



10 
11 
13 



*This document is a reformatted version of my PhD thesis. Advisor: Pro- 
fessor Victor M. Yakovenko. Committee: Professor J. Robert Dorfman, Pro- 
fessor Theodore L. Einstein, Professor Bei-Lok Hu, and Professor John D. 
Weeks. 



Acknowledgements 

My six year apprenticeship in physics at the University of 
Maryland was a unique and remarkable leg of my life. And to 
a large degree this is due to my advisor. Professor Victor M. 
Yakovenko, who has been at the center of my PhD education. 
I will deeply miss his clear perspective, his sure hand, and the 
incommensurate thrill of doing physics together. 



2 



I am grateful to Professors J. Robert Dorfman, Theodore 
L. Einstein, Bei-Lok Hu, and John D. Weeks for agreeing to 
serve on my advisory committee. 

Thanks are also due to Krishnendu Sengupta, Nicholas 
Dupuis, Hyok-Jon Kwon, Anatoley Zheleznyak, and Hsi- 
Sheng Goan all former students or post-docs of professor Vic- 
tor Yakovenko, for stimulating discussions in condensed mat- 
ter physics. My interests in physics were greatly molded by 
professor Bei-Lok Hu, through his remarkable courses in non- 
equilibrium physics, and also through personal interaction for 
more than three years. I had profited a lot from the association 
with Charis Anastopoulos and Sanjiv Shresta while working 
on quantum optics problems suggested by professor Bei-Lok 
Hu. 

It is not without melancholy that I thank my first true 
physics teachers: Ion Cotaescu, Ovidiu Lipan, Dan Luca, 
Adrian Neagu, Mircea Ra§a, Costel Rajinariu, Vlad Socol- 
iuc, and Dumitru Vulcanov. On an even deeper layer, I thank 
my family for all the support they provided me over the years, 
and for allowing me to pursue my interests. 

All the people mentioned above have generously shared 
with me some of their time and knowledge. I will do all I 
can to represent them well. 



I. STATISTICAL MECHANICS OF MONEY 
A. Introduction 

The application of statistical physics methods to economics 
promises fresh insights into problems traditionally not asso- 
ciated with physics (see, for example, the recent review and 
book'). Both statistical mechanics and economics study big 
ensembles: collections of atoms or economic agents, respec- 
tively. The fundamental law of equilibrium statistical mechan- 
ics is the Boltzmann-Gibbs law, which states that the proba- 
bility distribution of energy e is P{e) — Ce~^l^ , where T is 
the temperature, and C is a normalizing constant^. The main 
ingredient that is essential for the textbook derivation of the 
Boltzmann-Gibbs law^ is the conservation of energy^. Thus, 
one may generalize that any conserved quantity in a big sta- 
tistical system should have an exponential probability distri- 
bution in equilibrium. 

An example of such an unconventional Boltzmann-Gibbs 
law is the probability distribution of forces experienced by 
the beads in a cylinder pressed with an external force"*. Be- 
cause the system is at rest, the total force along the cylinder 
axis experienced by each layer of granules is constant and is 
randomly distributed among the individual beads. Thus the 
conditions are satisfied for the applicability of the Boltzmann- 
Gibbs law to the force, rather than energy, and it was indeed 
found experimentally"*. 

We claim that, in a closed economic system, the to- 
tal amount of money is conserved. Thus the equilibrium 
probability distribution of money P(m) should follow the 
Boltzmann-Gibbs law P(m) = Ce~™l^ . Here m is money, 
and T is an effective temperature equal to the average amount 
of money per economic agent. The conservation law of 



money^ reflects their fundamental property that, unlike mate- 
rial wealth, money (more precisely the fiat, "paper" money) is 
not allowed to be manufactured by regular economic agents, 
but can only be transferred between agents. Our approach 
here is very similar to that of Ispolatov et al.^. However, they 
considered only models with broken time-reversal symmetry, 
for which the Boltzmann-Gibbs law typically does not hold. 
The role of time-reversal symmetry and deviations from the 
Boltzmann-Gibbs law are discussed in detail in Sec. lIGI 

It is tempting to identify the money distribution P{m) with 
the distribution of wealth However, money is only one part 
of wealth, the other part being material wealth. Material prod- 
ucts have no conservation law because they can be manu- 
factured, destroyed, consumed, etc. Moreover, the monetary 
value of a material product (the price) is not constant. The 
same applies to stocks, which economics textbooks explicitly 
exclude from the definition of money^. So, we do not ex- 
pect the Boltzmann-Gibbs law for the distribution of wealth. 
Some authors believe that wealth is distributed according to a 
power law (Pareto-Zipf), which originates from a multiplica- 
tive random process^. Such a process may reflect, among 
other things, the fluctuations of prices needed to evaluate the 
monetary value of material wealth. 



B. Boltzmann-Gibbs distribution 

Let us consider a system of many economic agents 3> 1, 
which may be individuals or corporations. In this thesis, we 
only consider the case where their number is constant. Each 
agent i has some money and may exchange it with other 
agents. It is implied that money is used for some economic 
activity, such as buying or selling material products; however, 
we are not interested in that aspect. As in Ref.^, for us the 
only result of interaction between agents i and j is that some 
money A™ changes hands: [mi,mj] [m^,m'j] — [mi — 
Am, rrij + Am]. Notice that the total amount of money is 
conserved in each transaction: mi+mj ~ m'^+m'j. This local 
conservation law of money^ is analogous to the conservation 
of energy in coUisions between atoms. We assume that the 
economic system is closed, i. e. there is no external flux of 
money, thus the total amount of money M in the system is 
conserved. Also, in the first part of this chapter, we do not 
permit any debt, so each agent's money must be non-negative: 
TOi > 0. A similar condition applies to the kinetic energy of 
atoms: > 0. 

Let us introduce the probability distribution function of 
money P{m), which is defined so that the number of agents 
with money between m and m + dm is equal to NP{m) dm. 
We are interested in the stationary distribution P{m) corre- 
sponding to the state of thermodynamic equilibrium. In this 
state, an individual agent's money mi strongly fluctuates, but 
the overall probability distribution P{m) does not change. 

The equilibrium distribution function P{m) can be derived 
in the same manner as the equilibrium distribution function 
of energy P{e) in physics^. Let us divide the system into 
two subsystems 1 and 2. Taking into account that money 
is conserved and additive: m = mi + m2, whereas the 
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N=500, M=5*10^ time=4*10'^. 
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Money, m 

FIG. 1: Histogram and points: stationary probability distribution 
of money P(m). Solid curves: fits to the Boltzmann-Gibbs law 
P{m) (X exp(— m/T). Vertical lines: the initial distribution of 
money. 
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FIG. 2: Time evolution of entropy. Top curve: for the exchange of 
a random fraction u of the average money in the system: Am = 
1/ M/N . Bottom curve: for the exchange of a small constant amount 
Am = 1. The time scale for the bottom curve is 500 times greater 
than indicated, so it actually ends at the time 10^. 



C. Computer simulations 



probability is multiplicative: P = P1P2, we conclude that 
P{mi + 7112) = P{mi)P{m.2). The solution of this equa- 
tion is P{rn) — Ce^™l^\ thus the equilibrium probabil- 
ity distribution of money has the Boltzmann-Gibbs form. 
From the normalization conditions P{m) dm — 1 and 
mP{m)dm = M/N, we find that C = l/T and 
T = M/N . Thus, the effective temperature T is the aver- 
age amount of money per agent. 

The Boltzmann-Gibbs distribution can be also derived by 
maximizing the entropy of money distribution 

5 — — dmP{ni) InP(m) under the constraint of mon- 
ey conservation^ Following original Boltzmann's argument, 
let us divide the money axis < m < 00 into small bins of 
size dm and number the bins consecutively with the index b = 
1,2,... Let us denote the number of agents in a bin b as Nf,, 
the total number being N = J^bLi-^b- The agents in the bin 

6 have money mi,, and the total money is M — X^b^i "^f)-^6- 
The probability of realization of a certain set of occupation 
numbers { A*";,} is proportional to the number of ways N agents 
can be distributed among the bins preserving the set {Nt,}. 
This number is iV!/A^i!A^2! ■ • • The logarithm of probability 
is entropy IniV! — X^bli l^i^fc'- When the numbers Nf, are 
big and Stirling's formula In A^! w N\nN ~ N applies, the 
entropy per agent is 5 = (A^lniV - Y^bLi \nNb)/N = 
— Pb In Pfc, where Pt, = Nj, /N is the probability that an 
agent has money Wf,. Using the method of Lagrange multipli- 
ers to maximize the entropy S with respect to the occupation 
numbers {Ni,} with the constraints on the total money M and 
the total number of agents N generates the Boltzmann-Gibbs 
distribution for P{m)^. 



To check that these general arguments indeed work, we per- 
formed several computer simulations. Initially, all agents are 
given the same amount of money: P{m) = S{m — M/N), 
which is shown in Fig.[2as the double vertical line. One pair 
of agents at a time is chosen randomly, then one of the agents 
is randomly picked to be the "winner" (the other agent be- 
comes the "loser"), and the amount Am > is transferred 
from the loser to the winner. If the loser does not have enough 
money to pay {mi < Am), then the transaction does not take 
place, and we proceed to another pair of agents. Thus, the 
agents are not permitted to have negative money. This bound- 
ary condition is crucial in establishing the stationary distribu- 
tion. As the agents exchange money, the initial delta-function 
distribution first spreads symmetrically. Then, the probability 
density starts to pile up at the impenetrable boundary m — 0. 
The distribution becomes asymmetric (skewed) and ultimately 
reaches the stationary exponential shape shown in Fig.^ We 
used several trading rules in the simulations: the exchange of 
a small constant amount Am — 1, the exchange of a ran- 
dom fraction < < 1 of the average money of the pair: 
Am ~ v{mi + mj ) / 2, and the exchange of a random fraction 
V of the average money in the system: Am — v M/N . Fig- 
ures in this chapter mostly show simulations for the third rule; 
however, the final stationary distribution was found to be the 
same for all rules. 

In the process of evolution, the entropy S of the system 
increases in time and saturates at the maximal value for the 
BoltzmannGibbs distribution. This is illustrated by the top 
curve in Fig.|2]computed for the third rule of exchange. The 
bottom curve in Fig. |2] shows the time evolution of entropy 
for the first rule of exchange. The time scale for this curve is 
500 times greater than for the top curve, so the bottom curve 
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actually ends at the time 10^. The plot shows that, for the first 
rule of exchange, mixing is much slower than for the third one. 
Nevertheless, even for the first rule, the system also eventually 
reaches the Boltzmann-Gibbs state of maximal entropy, albeit 
over a time much longer than shown in Fig.|2] 

One might argue that the pairwise exchange of money may 
correspond to a medieval market, but not to a modern econ- 
omy. In order to make the model somewhat more realistic, 
we introduce firms. One agent at a time becomes a "firm". 
The firm borrows capital K from another agent and returns it 
with an interest rK, hires L agents and pays them wages W , 
manufactures Q items of a product and sells it to Q agents 
at a price R. All of these agents are randomly selected. The 
firm receives the profit F = RQ — LW — rK. The net re- 
sult is a many-body exchange of money that still satisfies the 
conservation law. 

Parameters of the model are selected following the pro- 
cedure described in economics textbooks^. The aggregate 
demand-supply curve for the product is taken to be R(ff) = 
V IQ^, where Q is the quantity people would buy at a price 
R, and r\ = 0.5 and V — 100 are constants. The produc- 
tion function of the firm has the conventional Cobb-Douglas 
form: Q{L, K) = Ll^K^^I^, where (3 = 0.8 is a constant. In 
our simulation, we set W = 10. By maximizing firm's profit 
F with respect to K and L, we find the values of the other 
parameters: L = 20, Q = 10, i? = 32, and F = 68. 

However, the actual values of the parameters do not matter. 
Our computer simulations show that the stationary probability 
distribution of money in this model always has the universal 
Boltzmann-Gibbs form independent of the model parameters. 



D. Thermal machine 

As explained in the Introduction, the money distribution 
P{m) should not be confused with the distribution of wealth. 
We believe that P{m) should be interpreted as the instanta- 
neous distribution of purchasing power in the system. Indeed, 
to make a purchase, one needs money. Material wealth nor- 
mally is not used directly for a purchase. It needs to be sold 
first to be converted into money. 

Let us consider an outside monopolistic vendor selling a 
product (say, cars) to the system of agents at a price p. Sup- 
pose that a certain small fraction / of the agents needs to buy 
the product at a given time, and each agent who has enough 
money to afford the price will buy one item. The fraction / 
is assumed to be sufficiently small, so that the purchase does 
not perturb the whole system significantly. At the same time, 
the absolute number of agents in this group is assumed to be 
big enough to make the group statistically representative and 
characterized by the Boltzmann-Gibbs distribution of money. 
The agents in this group continue to exchange money with the 
rest of the system, which acts as a thermal bath. The demand 
for the product is constantly renewed, because products pur- 
chased in the past expire after a certain time. In this situation, 
the vendor can sell the product persistently, thus creating a 
small steady leakage of money from the system to the vendor. 

What price p would maximize the vendor's income? 



To answer this question, it is convenient to introduce 
the cumulative distribution of purchasing power Af{m) — 
N P{m') dm' = iVe-"/^, which gives the number of 
agents whose money is greater than m. The vendor's income 
is fpJV{p). It is maximal when p = T,i. e. the optimal price 
is equal to the temperature of the system. This conclusion 
also follows from the simple dimensional argument that tem- 
perature is the only money scale in the problem. At the price 
p ~ T that maximizes the vendor's income, only the fraction 
J\f{T)/N =^ e"^ = 0.37 of the agents can afford to buy the 
product. 

Now let us consider two disconnected economic systems, 
one with the temperature Ti and another with T2: Ti > T2. 
A vendor can buy a product in the latter system at its equi- 
librium price T2 and sell it in the former system at the price 
Ti, thus extracting the speculative profit Ti — T2, as in a ther- 
mal machine. This example suggests that speculative profit 
is possible only when the system as a whole is out of equi- 
librium. As money is transferred from the high- to the low- 
temperature system, their temperatures become closer and 
eventually equal. After that, no speculative profit is possi- 
ble, which would correspond to the "thermal death" of the 
economy. This example brings to mind economic relations 
between developed and developing countries, with manufac- 
turing in the poor (low-temperature) countries for export to 
the rich (high-temperature) ones. 

We will demonstrate in Ch. HI] that for the large majority 
of the population the distribution of income is exponential. 
Hence, similar to the distribution of money, the distribution 
of income has a corresponding temperature. If in the previ- 
ous discussion about economic trade, the two countries have a 
different temperature for the distribution of income, the pos- 
sibility of constructing a thermal machine will still hold true. 
This is because purchasing power (total money) per unit time 
has two components: a positive one, income and a negative 
one, spending. Instead of making a one time purchase, the 
buyer will buy one product per unit time, effectively creating 
the thermal engine. We will revisit this idea in Sec. Ill El when 
we give the values for income temperature between the states 
of the USA and for the United Kingdom. 



E. Models with debt 

Now let us discuss what happens if the agents are permit- 
ted to go into debt. Debt can be viewed as negative money. 
Now when a loser does not have enough money to pay, he 
can borrow the required amount from a reservoir, and his 
balance becomes negative. The conservation law is not vi- 
olated: The sum of the winner's positive money and loser's 
negative money remains constant. When an agent with a 
negative balance receives money as a winner, she uses this 
money to repay the debt until her balance becomes positive. 
We assume for simplicity that the reservoir charges no inter- 
est for the lent money. However, because it is not sensible 
to permit unlimited debt, we put a limit nid on the maximal 
debt of an agent: to; > — m^. This new boundary condi- 
tion P{m < —md) = replaces the old boundary condition 
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FIG. 3: Histograms: stationary distributions of money with and with- 
out debt. Solid curves: fits to the Boltzmann-Gibbs laws with tem- 
peratures T = 1800 and T = f 000. 



P{m < 0) = 0. The result of a computer simulation with 
nid = 800 is shown in Fig.|3]together with the curve for = 
0. P{m) is again given by the Boltzmann-Gibbs law, but now 
with the higher temperature T — AI/N + nid, because the 
normalization conditions need to be maintained including the 
population with negative money: P{m) dm = 1 and 

m P{m) dm — M/N. The higher temperature makes 
the money distribution broader, which means that debt in- 
creases inequality between agents. 

In general, temperature is completely determined by the av- 
erage money per agent, (m) = M/N, and the boundary con- 
ditions. Suppose the agents are required to have no less money 
than mi and no more than m2- mi < m < m2- In this case, 
the two normalization conditions: f™^ P(m) dm = 1 and 

m P{m) dm = (m) with P{m) = C e^"'/"^ give the 
following equation for T 



tributions, the number of events becomes small and statistics 
poor, so the Boltzmann-Gibbs law loses applicability. Thus, 
we expect the Boltzmann-Gibbs law to hold only for the in- 
termediate range of money not too close either to the lower 
boundary or to the very high end. However, this range is the 
most relevant, because it covers the great majority of popula- 
tion. 

Lending creates equal amounts of positive (asset) and neg- 
ative (liability) money^. When economics textbooks de- 
scribe how "banks create money" or "debt creates money"^, 
they do not count the negative liabilities as money, and thus 
their money is not conserved. In our operational definition of 
money, we include all financial instruments with fixed denom- 
ination, such as currency, lOUs, and bonds, but not material 
wealth or stocks, and we count both assets and liabilities. With 
this definition, money is conserved, and we expect to see the 
Boltzmann-Gibbs distribution in equilibrium. Unfortunately, 
because this definition differs from economists' definitions of 
money (Ml, M2, M3, etc.^), it is not easy to find the appro- 
priate statistics. Of course, money can be also emitted by a 
central bank or government. This is analogous to an external 
influx of energy into a physical system. However, if this pro- 
cess is sufficiently slow, the economic system may be able to 
maintain quasi-equilibrium, characterized by a slowly chang- 
ing temperature. 

We performed a simulation of a model with one bank and 
many agents. The agents keep their money in accounts on 
which the bank pays interest. The agents may borrow money 
from the bank, for which they must pay interest in monthly 
installments. If they cannot make the required payments, they 
may be declared bankrupt, which relieves them from the debt, 
but the liability is transferred to the bank. In this way, the con- 
servation of money is maintained. The model is too elaborate 
to describe it in full detail here. We found that, depending on 
the parameters of the model, either the agents constantly lose 
money to the bank, which steadily reduces the agents' temper- 
ature, or the bank constantly loses money, which drives down 
its own negative balance and steadily increases the agents' 
temperature. 



coth 



Am 



T 
Am 



(1) 



wherem= (mi+m2)/2 and Am = (m2— mi)/2. It follows 
from Eq. that the temperature is positive when m > (m), 
negative when m < (m), and infinite (P(m) ~ const) when 
m — (m). In particular, if agents' money are bounded from 
above, but not from below: — oo < m < m2, the temper- 
ature is negative. That means an inverted Boltzmann-Gibbs 
distribution with more rich agents than poor. 

Imposing a sharp cutoff at m^ may be not quite realistic. 
In practice, the cutoff may be extended over some range de- 
pending on the exact bankruptcy rules. Over this range, the 
Boltzmann-Gibbs distribution would be smeared out. So we 
expect to see the Boltzmann-Gibbs law only sufficiently far 
from the cutoff region. Similarly, in experimenf*, some devi- 
ations from the exponential law were observed near the lower 
boundary of the distribution. Also, at the high end of the dis- 



F. Boltzmann equation 

The Boltzmann-Gibbs distribution can be also derived from 
the Boltzmann equation^, which describes the time evolution 
of the distribution function P{m) due to pairwise interactions: 



^]P(m)P(m') (2) 



dP{m) 
Jt ^ 

+W[„i_A,m'+A]^[m,m']-P(n^ - A)P{m' + A)} dm' dA. 

Here W[,n,m']^[,n-A,m'+A] is the rate of transferring money 
A from an agent with money m to an agent with money m'. 
If a model has time-reversal symmetry, then the transition rate 
of a direct process is the same as the transition rate of the re- 
versed process, thus the ui-factors in the first and second lines 
of Eq. (|3 are equal. In this case, it can be easily checked that 
the Boltzmann-Gibbs distribution P{m) = C exp{—m/T) 
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N=500, M=5*10^ a=1/3. 
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FIG. 4: Histogram: stationary probability distribution of money in 
the multiplicative random exchange model studied in Ref.*. Solid 
curve: the Boltzmann-Gibbs law. 



N=500, IV1=5*10^ tax=40%. 

16| 1 1 1 1 — 




500 1000 1500 2000 2500 3000 

Money, m 



FIG. 5: Histogram: stationary probability distribution of money in 
the model with taxes and subsidies. Solid curve: the Boltzmann- 
Gibbs law. 



nullifies the right-hand side of Eq. (|2j; thus this distribution 
is stationary: dP{m)/dt — G-. 



G. Non-Boltzmann-Gibbs distributions 

However, if time-reversal symmetry is broken, the two tran- 
sition rates w in Eq. (|3 may be different, and the system may 
have a non-Boltzmann-Gibbs stationary distribution or no sta- 
tionary distribution at all. Examples of this kind were stud- 
ied in Ref.^. One model was called multiplicative random 
exchange. In this model, a randomly selected loser i loses 
a fixed fraction a of his money rrii to a randomly selected 
winner j: [mi,mj] [(1 — a)mi , mj + ami]. If we try 
to reverse this process and immediately appoint the winner 
i to become a loser, the system does not return to the orig- 
inal configuration [mi,mj\: [(1 — a)mi , rrij + ami] 
[(1 — a)mi + a(mj + ami) , (1 — a)(mj + ami)]. Ex- 
cept for a — 1/2, the exponential distribution function is not 
a stationary solution of the Boltzmann equation derived for 
this model in Ref.''. Instead, the stationary distribution has the 
shape shown in Fig.|4]for a = 1/3, which we reproduced in 
our numerical simulations. It still has an exponential tail end 
at the high end, but drops to zero at the low end for a < 1/2. 
Another example of similar kind was studied in Refi", which 
appeared after the first version of our paper was posted as 
rcond-mat/0001432 on January 30, 2000. In that model, the 
agents save a fraction A of their money and exchange a ran- 
dom fraction e of their total remaining money: [mi, mj] 
[Ami + e(l — X){mi + mj) , Xmj + (1 — e)(l — \){mi + mj)]. 
This exchange also does not return to the original configura- 
tion after being reversed. The stationary probability distribu- 
tion was found in Ref.'" to be nonexponential for A 7^ with 
a shape qualitatively similar to the one shown in Fig.O 

Another interesting example which has a non-Boltzmann- 



Gibbs distribution occurs in a model with taxes and subsi- 
dies. Suppose a special agent ("government") collects a frac- 
tion ("tax") of every transaction in the system. The collected 
money is then equally divided among all agents of the sys- 
tem, so that each agent receives the subsidy Sm with the fre- 
quency 1/ts. Assuming that Sm is small and approximating 
the collision integral with a relaxation time we obtain the 
following Boltzmann equation 

dP{m) ^ 6m dP{m) _ P{m) - P{m) 
dt Ts dm Tr ' 

where P{m) is the equilibrium Boltzmann-Gibbs function. 
The second term in the left-hand side of Eq. Q is analo- 
gous to the force applied to electrons in a metal by an ex- 
ternal electric field^. The approximate stationary solution 
of Eq. Q is the displaced Boltzmann-Gibbs one: P{m) = 
P{m—{Tr/Ts) Sm). The displacement of the equilibrium dis- 
tribution P{m) by (Tr/Ts) 6m would leave an empty gap near 
TO = 0. This gap is filled by interpolating between zero pop- 
ulation at m = and the displaced distribution. The curve 
obtained in a computer simulation of this model (Fig.|3 qual- 
itatively agrees with this expectation. The low-money popula- 
tion is suppressed, because the government, acting as an exter- 
nal force, "pumps out" that population and pushes the system 
out of thermodynamic equilibrium. We found that the entropy 
of the stationary state in the model with taxes and subsidies is 
few a percent lower than without. 

These examples show that the Boltzmann-Gibbs distribu- 
tion is not fully universal, meaning that it does not hold for 
just any model of exchange that conserves money. Neverthe- 
less, it is universal in a limited sense: For a broad class of 
models that have time-reversal symmetry, the stationary dis- 
tribution is exponential and does not depend on the details of 
a model. Conversely, when time-reversal symmetry is bro- 
ken, the distribution may depend on model details. The dif- 
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ference between these two classes of models may be rather 
subtle. For example, let us change the multiplicative random 
exchange from a fixed fraction of loser's money to a fixed 
fraction of the total money of winner and loser. This mod- 
ification retains the multiplicative idea that the amount ex- 
changed is proportional to the amount involved, but restores 
time-reversal symmetry and the Boltzmann-Gibbs distribu- 
tion. In the model with Am = 1 discussed in the next Section, 
the difference between time-reversible and time-irreversible 
formulations amounts to the difference between impenetra- 
ble and absorbing boundary conditions at to = 0. Unlike in 
physics, in economy there is no fundamental requirement that 
interactions have time-reversal symmetry. However, in the ab- 
sence of detailed knowledge of real microscopic dynamics of 
economic exchange, the semiuniversal Boltzmann-Gibbs dis- 
tribution appears to be a natural starting point. 

Moreover, deviations from the Boltzmann-Gibbs law may 
occur only if the transition rates w in Eq. (|3 explicitly de- 
pend on the agents' money m or m' in an asymmetric man- 
ner. In another simulation, we randomly preselected winners 
and losers for every pair of agents In this case, money 

flows along directed links between the agents: i ^ j ^ k, 
and time-reversal symmetry is strongly broken. This model 
is closer to the real economy, in which, for example, one 
typically receives money from an employer and pays it to a 
grocer, but rarely the reverse. Nevertheless, we still found 
the Boltzmann-Gibbs distribution of money in this model, be- 
cause the transition rates w do not explicitly depend on m and 
to'. 



H. Nonlinear Boltzmann equation vs. linear master equation 

For the model where agents randomly exchange the con- 
stant amount Am = 1, the Boltzmann equation is: 

CO oo 

— Pm+l Pn + Pm-1 Pn 

n—0 n—1 

oo oc 

- Prn Y.^-- ^™ (4) 

n—0 n—1 
= (-Pm+l + -Pr?j-1 ^ 2P„j) + Pq{P„i — Pm-i), (5) 

where Pm = P{m) and we have used X]m=o ~ 1- 
first, diffusion term in Eq. (|5} is responsible for broadening of 
the initial delta-function distribution. The second term, pro- 
portional to Pq^ is essential for the Boltzmann-Gibbs distri- 
bution P,„ = e~™/-^(l — e^^/-^) to be a stationary solution 
of Eq. (|5}- In a similar model studied in Ref.'', the second 
term was omitted on the assumption that agents who lost all 
money are eUminated: Pq = 0. In that case, the total number 
of agents is not conserved, and the system never reaches any 
stationary distribution. Time-reversal symmetry is violated, 
since transitions into the state to = are permitted, but not 
out of this state. 

If we treat Pq a constant, Eq. (|5jl looks like a linear 
Fokker- Planck equation"^ for P^, with the first term describing 
diffusion and the second term an external force proportional 



to Pq. Similar equations were studied in Ref.**. Eq. Q can be 
also rewritten as 

dP 

= P,„+i - (2 - Po)P,„ + (1 - Po)P™-i. (6) 

The coefficient ( 1 — Pq) in front of Pm - 1 represents the rate of 
increasing money by Am — 1, and the coefficient 1 in front of 
Pm+i represents the rate of decreasing money by Am = — 1. 
Since Po > 0, the former is smaller than the latter, which 
results in the stationary Boltzmann-Gibbs distributions P„j ~ 
(1 — Pq)™. An equation similar to Eq. (|6j describes a Markov 
chain studied for strategic market games in Ref. ' ' . Naturally, 
the stationary probability distribution of wealth in that model 
was found to be exponential". 

Even though Eqs. Q and (|6j look like linear equations, 
nevertheless the Boltzmann equation (|2j and is a pro- 
foundly nonlinear equation. It contains the product of two 
probability distribution functions P in the right-hand side, be- 
cause two agents are involved in money exchange. Most stud- 
ies of wealth distribution^ have the fundamental flaw that they 
use a single-particle approach. They assume that the wealth 
of an agent may change just by itself and write a linear mas- 
ter equation for the probability distribution. Because only 
one particle is considered, this approach cannot adequately 
incorporate conservation of money. In reality, an agent can 
change money only by interacting with another agent, thus the 
problem requires a two-particle probability distribution func- 
tion. Using Boltzmann's molecular chaos hypothesis, the two- 
particle function is factorized into a product of two single- 
particle distributions functions, which results in the nonlinear 
Boltzmann equation. Conservation of money is adequately in- 
corporated in this two-particle approach, and the universality 
of the exponential Boltzmann-Gibbs distribution is transpar- 
ent. 

I. Conclusions 

Everywhere in this chapter we assumed some randomness 
in the exchange of money. Our results would apply the best to 
the probability distribution of money in a closed community 
of gamblers. In more traditional economic studies, the agents 
exchange money not randomly, but following deterministic 
strategies, such as maximization of utility functions^- The 
concept of equilibrium in these studies is similar to mechan- 
ical equilibrium in physics, which is achieved by minimizing 
energy or maximizing utility. However, for big ensembles, 
statistical equilibrium is a more relevant concept. When many 
heterogeneous agents deterministically interact and spend var- 
ious amounts of money from very little to very big, the money 
exchange is effectively random. In the future, we would like 
to uncover the Boltzmann-Gibbs distribution of money in a 
simulation of a big ensemble of economic agents following re- 
alistic deterministic strategies with money conservation taken 
into account. That would be the economics analog of molecu- 
lar dynamics simulations in physics. While atoms collide fol- 
lowing fully deterministic equations of motion, their energy 
exchange is effectively random due to the complexity of the 
system and results in the Boltzmann-Gibbs law. 
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We do not claim that the real economy is in equilibrium. 
(Most of the physical world around us is not in true equilib- 
rium either.) Nevertheless, the concept of statistical equilib- 
rium is a very useful reference point for studying nonequilib- 
rium phenomena. 



n. DISTRIBUTION OF INCOME AND WEALTH 

A. Introduction 

In Ch.|I]we predicted that the distribution of money should 
follow an exponential Boltzmann-Gibbs law. Unfortunately, 
we were not able to find data on the distribution of money. On 
the other hand, we found many sources with data on income 
distribution for the United States (USA) and United Kingdom 
(UK), as well as data on the wealth distribution in the UK, 
which are presented in this chapter. In all of these data, we 
find that the great majority of the population is described by 
an exponential distribution, whereas the high-end tail follows 
a power law. 

The study of income distribution has a long history. 
Paretc''* proposed in 1897 that income distribution obeys a 
universal power law valid for all times and countries. Sub- 
sequent studies have often disputed this conjecture. In 1935, 
Shirras- concluded: "There is indeed no Pareto Law. It is 
time it should be entirely discarded in studies on distribu- 
tion". Mandelbrot'^ proposed a "weak Pareto law" applica- 
ble only asymptotically to the high incomes. In such a form, 
Pareto's proposal is useless for describing the great majority 
of the population. 

Many other distributions of income were proposed: Levy, 
log-normal, Champernowne, Gamma, and two other forms by 
Pareto himself (see a systematic survey in the World Bank re- 
search publication'^). Theoretical justifications for these pro- 
posals form two schools: socio-economic and statistical. The 
former appeals to economic, political, and demographic fac- 
tors to explain the distribution of income (e. g."*), whereas the 
latter invokes stochastic processes. Gibrat'^ proposed in 1931 
that income is governed by a multiplicative random process, 
which results in a log-normal distribution (see alsc^"). How- 
ever, Kalecki^' pointed out that the width of this distribution 
is not stationary, but increases in time. Levy and Solomon^^ 
proposed a cut-off at lower incomes, which stabilizes the dis- 
tribution to a power law. 

Many researchers tried to deduce the Pareto law from a the- 
ory of stochastic processes. Gibrat'^ proposed in 1931 that 
income and wealth are governed by multiplicative random 
processes, which result in a log-normal distribution. These 
ideas were later followed, among many others, by Montroll 
and Shlesingej2£. However, already in 1945 Kalecki^ pointed 
out that the log-normal distribution is not stationary, because 
its width increases in time. Modern econophysicists^^-^^-^"' 
also use various versions of multiplicative random processes 
in theoretical modeling of wealth and income distributions. 

Unfortunately, numerous recent papers on this subject do 
very little or no comparison at all with real statistical data, 
much of which is widely available these days on the Inter- 



net. In order to fill this gap, we analyzed the data on in- 
come distribution in the United States (US) from the Bureau 
of Census and the Internal Revenue Service (IRS) in Refi^. 
We found that the individual income of about 95% of popu- 
lation is described by the exponential law. The exponential 
law, also known in physics as the Boltzmann-Gibbs distribu- 
tion, is characteristic for a conserved variable, such as energy. 
In Ref.^'', we argued that, because money (cash) is conserved, 
the probability distribution of money should be exponential. 
Wealth can increase or decrease by itself, but money can only 
be transferred from one agent to another. So, wealth is not 
conserved, whereas money is. The difference is the same as 
the difference between unrealized and realized capital gains in 
stock market. 

In Sec. Ill Bl we propose that the distribution of individual 
income is given by an exponential function. This conjecture 
is inspired by the results of Ch.lH^ where we argued that the 
probability distribution of money in a closed system of agents 
is given by the exponential Boltzmann-Gibbs function. We 
compare our proposal with the census and tax data for individ- 
ual income in USA. In Sec. Ill Dl we show that the exponential 
distribution has to be amended for the top 5% of incomes by 
a power-law function. In Sec. Ill Fl we derive the distribution 
function of income for families with two earners and compare 
it with census data. In Sec. Ill Gl we discuss the distribution of 
wealth. In Sec. Ill HI we critically examine several other alter- 
natives proposed in the literature for the distribution of income 
and wealth. Speculations on the possible origins of the expo- 
nential and power-law distribution of income and conclusions 
for this chapter are given in Sec. Ill II 



B. Distribution of income for individuals 

C. Exponential distribution of income 

We denote income by the letter r (for "revenue"). The prob- 
ability distribution function of income, P(r), (called the prob- 
ability density in book'^) is defined so that the fraction of in- 
dividuals with income between r and r + dr is P{r) dr. This 
function is normalized to unity (100%): P[f) dr = 1. We 
propose that the probability distribution of individual income 
is exponential: 

Pi(r) =cxp(-r/i?)/i?, (7) 

where the subscript 1 indicates individuals. Function 
contains one parameter R, equal to the average income: 
/(J r Pi (r) dr = B, and analogous to temperature in the 
Boltzmann-Gibbs distribution^^. 

From the Survey of Income and Program Participation 
(SIPP)2^ we downloaded the variable TPTOINC (total in- 
come of a person for a month) for the first "wave" (a four- 
month period) in 1996. Then we eliminated the entries with 
zero income, grouped the remaining entries into bins of the 
size 10/3 k$, counted the numbers of entries inside each bin, 
and normalized to the total number of entries. The results are 
shown as the histogram in Fig.|6j where the horizontal scale 
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Individual annual income, k$ 

FIG. 6: Histogram: Probability distribution of individual income 
from the U.S. Census data for 1996"^. Solid line: Fit to the expo- 
nential law. Inset plot A: The same with the logarithmic vertical 
scale. Inset plot B: Cumulative probability distribution of individual 
income from PSID for 199321. 

has been multiplied by 12 to convert monthly income to an 
annual figure. The sohd line represents a fit to the exponential 
function 0. In the inset, plot A shows the same data with 
the logarithmic vertical scale. The data fall onto a straight 
line, whose slope gives the parameter R in Eq. Q. The ex- 
ponential law is also often written with the bases 2 and 10: 
Pi(r) (X 2-''/^=^ (X lO"''/-^". The parameters R, R2 and 
Rio are given in line (c) of TableU 





Source 


Year 


R($) 


R2 ($) 


-R10 ($) 


Set size 


a 


PSIDii 


1992 


18,844 


13,062 


43,390 


1.39x10^ 


b 


IRS^ 


1993 


19,686 


13,645 


45,329 


1.15x10** 


c 


SIPPjjS 


1996 


20,286 


14,061 


46,710 


2.57x10^ 


d 


SIPPfS 


1996 


23,242 


16,110 


53,517 


1.64x10^ 


e 


IRS^' 


1997 


35,200 


24,399 


81,051 


1.22x10® 



TABLE I: Parameters R, R2, and Rio obtained by fitting data from 
different sources to the exponential law with the bases e, 2, and 
10, and the sizes of the statistical data sets. 

Plot B in the inset of Fig. |6l shows the data from the Panel 
Study of Income Dynamics (PSID) conducted by the Insti- 
tute for Social Research of the University of Michigan^. We 
downloaded the variable V30821 "Total 1992 labor income" 
for individuals from the Final Release 1993 and processed the 
data in a similar manner. Shown is the cumulative probabil- 
ity distribution of income N{r) (called the probability dis- 
tribution in book'^). It is defined as N{r) = J^°° P[r') dr' 
and gives the fraction of individuals with income greater than 
r. For the exponential distribution Q, the cumulative dis- 
tribution is also exponential: Ni{r) — Pi{r')dr' = 
exp(— r/i?). Thus, R2 is the median income; 10% of popula- 
tion have income greater than i?io and only 1 % greater than 
2i?io. The points in the inset fall onto a straight line in the 
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FIG. 7: Points: Cumulative fraction of tax returns vs income from 
the IRS data for 1997"'. Solid line: Fit to the exponential law. Inset 
plot A: The same with the logarithmic vertical scale. Inset plot B: 
Probability distribution of individual income from the IRS data for 
I9932&. 

logarithmic scale. The slope is given in line (a) of TableU 

The points in Fig. show the cumulative distribution of 
tax returns vs income in 1997 from column 1 of Table 1.1 
of Ref.^^. (We merged 1 k$ bins into 5 k$ bins in the interval 
1-20 k$.) The solid line is a fit to the exponential law. Plot A 
in the inset of Fig.0shows the same data with the logarithmic 
vertical scale. The slope is given in line (e) of Table U Plot 
B in the inset of Fig. shows the distribution of individual 
income from tax returns in 1993^^. The logarithmic slope is 
given in line (b) of Table|l] 

While Figs. |5] and clearly demostrate the fit of income 
distribution to the exponential form, they have the following 
drawback. Their horizontal axes extend to +00, so the high- 
income data are left outside of the plots. The standard way to 
represent the full range of data is the so-called Lorenz curve 
(for an introduction to the Lorenz curve and Gini coefficient, 
see book'^). The horizontal axis of the Lorenz curve, x{r), 
represents the cumulative fraction of population with income 
below r, and the vertical axis y{r) represents the fraction of 
income this population accounts for: 

xir) = / F(r )dr , v(r) — — ^ . (8) 

^ ^ Jo J^r'P{r')dr' 

As r changes from to 00, x and y change from to 1, and 
Eq. (jS) parametrically defines a curve in the (x, j/)-space. 
Substituting Eq. into Eq. (|8}, we find 

x{f) = 1 — exp(— f), y{f) — x{f) — rexp(— 7~), (9) 

where f = r/ R. Excluding f, we find the explicit form of the 
Lorenz curve for the exponential distribution: 

y^x+{l-x)\n{l-x). (10) 

R drops out, so Eq. dlOt has no fitting parameters. 
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FIG. 8: Solid curve: Lorenz plot for the exponential distribution. 
Points: IRS data for 1979-199'?2Z. Inset points: Gini coefficient data 
from IRS^'. Inset line: The calculated value 1/2 of the Gini coeffi- 
cient for the exponential distribution. 

The function jlO> is shown as the solid curve in Fig.|S] The 
straight diagonal line represents the Lorenz curve in the case 
where all population has equal income. Inequality of income 
distribution is measured by the Gini coefficient G, the ratio of 
the area between the diagonal and the Lorenz curve to the area 
of the triangle beneath the diagonal 

G = 2 [ {x-y)dx (11) 

The Gini coefficient is confined between (no inequality) and 
1 (extreme inequality). By substituting Eq. ilO\ into the inte- 
gral, we find the Gini coefficient for the exponential distribu- 
tion: Gi = 1/2. 

The points in Fig. |8] represent the tax data during 1979- 
1997 fromRef.^^. With the progress of time, the Lorenz points 
shifted downward and the Gini coefficient increased from 0.47 
to 0.56, which indicates increasing inequality during this pe- 
riod. However, overall the Gini coefficient is close to the value 
0.5 calculated for the exponential distribution, as shown in the 
inset of Fig.|8l 

D. Power-law tail and "Bose" condensation 

As Fig.lSjshows, the Lorenz curve deviates from the theo- 
retical Lorenz curve implied by the exponential distribution, 
mostly for the top 20% of tax returns. Moreover, as explained 
in the previous section, the Lorenz curve for a pure exponen- 
tial distribution is independent of temperature (the scale of the 
distribution). Therefore the variations in the Lorenz curves 
over the period 1979-1997 suggest that the shape of the distri- 
bution changes. 
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FIG. 9: Cumulative probability distribution of US individual income 
for 1997 in log-log scale, with points (raw data) and solid lines (ex- 
ponential and power-law fit). The inset shows the exponential regime 
and the fit with a Boltzmann-Gibbs distribution in the log-linear 
scale. 

In Sec. lll Al we mentioned the early proposal of Pareto, who 
claimed that the income distribution obeys a universal power 
law valid for all times and countriesi^. The Pareto distribution 
P{r) — A/r" has several undesirable properties. It diverges 
for low incomes, and if a < 2 the distribution has a diverg- 
ing first moment. These two properties are never found in real 
data. Moreover, our study has shown that at least for the over- 
whelming majority (95%) of the population the distribution of 
income is exponential. 

The data for the remaining top 5% percent of the population 
is hard to get and often it is unrehable. Because the number 
of people with income in the top few percents is small, census 
data is poor. The Internal Revenue Service (IRS) has reliable 
data but it seldom reports it. We managed to find a data set 
from IRS which contains high income data^^. 

Fig.|9]shows United States IRS income data for 19972^. The 
data points represent the cumulative probability distribution in 
log-log scale. Although for large income values r > 120k$ 
we have only three data points, the points align on a straight 
line in accordance with a Pareto power-law distribution. 

From the plot it is evident that there is a discontinuity 
around 120k$ where the exponential and power-law regimes 
intersect. We conclude that it is not possible to describe 
the entire income distribution with one single differentiable 
function. The two regimes of the probability distribution 
are clearly separated and this may be due to different in- 
come dynamics in the two regimes. Among others, the 
Adjusted Gross Income contains capital gains that is stock- 
market gains/losses. It is conceivable that for the top 5% of 
the population, capital gains rather than the labor income ac- 
count for the majority share of the total income. The capital 
gains contribution to the Adjusted Growth Income may be re- 
sponsible for the observed power-law in the tail of the income 
distribution. 

We plot the Lorenz curve for income in Fig.^| As in Sec. 
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FIG. 10: Lorenz plot with points (raw data) and solid line (function 
in\ witii fraction b = 16%). 

Ill B I the horizontal and vertical coordinates are the cumula- 
tive population x{r) and the cumulative income y{r) from (|8|i. 
An imaginary line through the data shows an abrupt change 
in derivative for the last 2% of the tax returns. This sudden 
change in derivative around the 98% mark on the horizontal 
axis of Fig. ^1 gives an almost infinite slope for the Lorenz 
curve. Physically, an infinite slope in the Lorenz function for 
arguments x just slightly less than one, is equivalent to hav- 
ing a finite amount of total society's income in the hand of a 
very few individuals. We have coined for this effect the name 
"Bose condensation of income" because of similarities with 
the celebrated phase-transition from statistical mechanics. 

To understand quantitatively the "Bose condensation" of 
income, we calculate / the ratio of the total income of the 
society if all the population is described by the exponential 
law, to the actual total income of the population. For the 
USA tax income data for 1997 shown in Fig. ^| this frac- 
tion was / ~ 0.84, which means that the "condensate" has 
6=1—/ = 16% of the total income of the population. The 
analytical formula for the Lorenz curve in this case is 

y= {l-b)[x + {l-x)lii{l-x)]+bdil~x), (12) 

where the first term represents the contribution of the 
Boltzmann-Gibbs exponential regime and the second term 
represents the contribution of the Pareto power-law tail. It is 
remarkable that as in the case of ilO\ . Ea.ll2ldoes not have any 
fitting parameters. The condensate fraction b is completely de- 
termined by temperature, which in turn is determined from the 
probability distribution. 

The function il2i is plotted as a solid line in Fig.|5] One can 
see that the data systematically deviate from the exponential 
law because of the income concentrated in the power-law tail; 
however, the deviation is not very big. The inequality of the 
US income distribution was also increasing during that time 
period^^, which implies that the size of the "Bose condensate" 
has increased in time, too. 
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FIG. 11: Cumulative probability distributions of yearly individual 
income for different states of the USA shown as raw data (top inset) 
and scaled data in log-log, log-linear (lower inset). 

The Gini coefficient defined in ( II 1> can be calculated for the 
Lorenz curve given by (I12> . The effect of the power-law tails 
changes the value of the Gini coefficient from Gi = 1/2 in 
the case of a pure exponential distribution to Gf, = (1 + b)/2. 



E. Geographical variations in income distribution 

In the previous sections we found that the distribution of 
income is an exponential for the large majority of the popu- 
lation followed by a power-law tail for the top few percent of 
the population. It would be interesting to estabhsh the univer- 
sality of this shape for various other countries. 

We obtained the data on distribution of the yearly individ- 
ual income in 1998 for each of the 50 states and the District 
of Columbia that constitute the USA from the Web site of the 
IRSi^. We plot the original raw data for the cumulative distri- 
bution of income in log-linear scale in the upper inset of Fig. 
rm The points spread significantly, particularly at higher in- 
comes. For example, the fraction of individuals with income 
greater than 1 M$ varies by an order of magnitude between 
different states. However, after we rescale the data in the man- 
ner described in the preceding section, the points collapse on 
a single curve shown in log-log scale in the main panel and 
log-linear scale in the lower inset. The open circles represent 
the US average, obtained by treating the combined data for 
all states as a single set. We observe that the distribution of 
higher incomes approximately follows a power law with the 
exponent a = 1.7 ± 0.1, where the ±0.1 variation includes 
70% of all states. On the other hand, for about 95% of indi- 
viduals with lower incomes, the distribution follows an expo- 
nential law with the average US temperature Rjjs = 36.4 k$ 
(r'i^I = 25.3 k$ and = 83.9 k$). The temperatures of 

the individual states differ from Rjjs by the amounts shown in 
Tableim For example, the temperature of Connecticut (CT) is 
1 .25 times higher and the temperature of West Virginia (WV) 




1 00% 



2 1 0% 

X 

to 



1% 



0.1% 



100%C>- 



10% 




*4l 

I 40 80 120 

Adjusted gross Income, k$ 



12 



All states of the USA, IRS data for 1 998 
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FIG. 12: Lorenz plot with points (raw data) and solid line (function 
<10> calculated for the exponential law). 



TABLE II: Deviations of the state temperatures from the average US 
temperature. 



CT 


NJ 


MA 


MD 


VA 


CA 


NY 


IL 


CO 


NH 


AK 


25% 


24% 


14% 


14% 


9% 


9% 


7% 


6% 


6% 


5% 


5% 


DC 


DE 


Ml 


WA 


MN 


GA 


TX 


Rl 


AZ 


PA 


FL 


5% 


4% 


4% 


2% 


1% 


0% 


-1% 


-3% 


-3% 


-3% 


-4% 


KS 


OR 


HI 


NV 


NC 


WI 


IN 


UT 


MO 


VT 


TN 


-5% 


-6% 


-7% 


-7% 


-7% 


-8% 


-8% 


-9% 


-9% 


-9% 


-11% 


NE 


OH 


LA 


AL 


SC 


lA 


WY 


NM 


KY 


ID 


OK 


-12% 


-12% 


-13% 


-13% 


-13% 


-14% 


-14% 


-14% 


-14% 


-15% 


-16% 


ME 


MT 


AR 


SD 


ND 


MS 


WV 










-16% 


-19% 


-19% 


-20% 


-20% 


-21% 


-22% 











is 0.78 times lower than the average US temperature. 

The Lorenz plot for all states is shown in Fig. 1121 together 
with the solid curve representing Eq. ilQ\ . The majority of 
points are well clustered and are not too far from the solid 
curve. The exceptions are Wyoming (WY) with much higher 
inequality and the Washington state (WA) with noticeably 
lower inequality of income distribution. The average US 
data, shown by open circles, is consistent with our previous 
results^^. Unlike in the UK case, we did not make any ad- 
justment in the US case for individuals with income below the 
threshold, which appears to be sufficiently low. 

We obtained the data on the yearly income distribution in 
the UK for 1997/8 and 1998/9 from the Web site of the IR^^. 
The data for 1994/5, 1995/6, and 1996/7 were taken from the 
Annual Abstract of Statistics derived from the IR"*". The data 
for these 5 years are presented graphically in Fig.^] In the 
upper inset, the original raw data for the cumulative distribu- 
tion are plotted in log-linear scale. For not too high incomes, 
the points form straight lines, which implies the exponential 
distribution N{r) cx exp(— r/i?), where r stands for income 
(revenue), and R is the income "temperature". However, the 
slopes of these lines are different for different years. The tem- 
peratures for the years 1994/5, 1995/6, 1996/7, and 1997/8 
differ from the temperature for 1998/9, = 11.7 ki: 
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FIG. 13: Cumulative probability distributions of yearly individual 
income in the UK shown as raw data (top inset) and scaled data in 
log-log (left panel), log-linear (lower inset), and Lorenz (right panel) 
coordinates. Solid curve: fit to function jlOt calculated for the expo- 
nential law. 



(R^Jl^ = 8.1 k£ and R'j^"^^ = 26.9 k£), by the factors 0.903, 
0.935, 0.954, and 0.943. To compensate for this effect, we 
rescale the data. We divide the horizontal coordinates (in- 
comes) of the data sets for different years by the quoted above 
factors and plot the results in log-log scale in the main panel 
and log-linear scale in the lower inset. We observe scaling: the 
collapse of points onto a single curve. Thus, the distributions 
Ni (r) for different years i are described by a single function 
f{r/Ri). The main panel shows that this scaling function / 
follows a power law with the exponent a = 2.0-2.3 at high in- 
comes. The lower inset shows that / has an exponential form 
for about 95% of individuals with lower incomes. These re- 
sults qualitatively agree with a similar study by Cranshaw'^i. 
He proposed that the P{r) data for lower incomes are better 
fitted by the Gamma distribution T[r) oc exp(— r/i?). For 
simplicity, we chose not to introduce the additional fitting pa- 
rameter p. 

We must mention that the individuals with income below a 
certain threshold are not required to report to the IR. That is 
why the data in the lower inset do not extend to zero income. 
We extrapolate the straight line to zero income and take the 
intercept with the vertical axis as 100% of individuals. Thus, 
we imply that the IR data does not account for about 25-27% 
of individuals with income below the threshold. 

The Lorenz curve for the distribution of the UK income 
is shown in Fig.^J We treat the number of individuals be- 
low the income threshold and their total income as adjustable 
parameters, which are the horizontal and vertical offsets of 
the coordinates origin relative to the lowest known data point. 
These parameters are chosen to fit the Lorenz curve for the 
exponential law dlOt shown as the solid line. The fit is very 
good, and is illustrated in Fig.^] The horizontal offsets are 
28-34%, which is roughly consistent with the numbers quoted 
for the lower inset of the left panel. 
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FIG. 14: Cumulative probability distributions of yearly individual 
income in the UK shown as raw data (top inset) and scaled data in 
log-log (left panel), log-linear (lower inset), and Lorenz (right panel) 
coordinates. Solid curve: fit to function <10> calculated for the expo- 
nential law. 

The income temperature for the UK in 1998/9 was Ruk = 
11.7 k£ and for the US in 1998 was Rus = 36.4 k$. Us- 
ing the exchange rate as of December 31, 1998 to convert 
pounds into dollars'*^, we find that the UK temperature was 
Ruk — 19.5 k$, which is 1.87 times lower than the US tem- 
perature. The difference in temperatures indicates nonequilib- 
rium, which can be exploited to create a thermal machine^*. 
The gain (profit) produced by such a thermal machine is pro- 
portional to the difference in temperatures. In agreement with 
the second law of thermodynamics, money would flow from a 
high-temperature system to a low-temperature one. This may 
explain the huge trade deficit of the USA in global interna- 
tional trade with other, lower-temperature countries. The vari- 
ation of temperatures between different states of the USA is 
shown in Table HII 

F. Distribution of income for 
families 

Now let us discuss the distribution of income for families 
with two earners. The family income r is the sum of two indi- 
vidual incomes: r = ri+r2- Thus, the probability distribution 
of the family income is given by the convolution of the indi- 
vidual probability distributions^^. If the latter are given by the 
exponential function Q, the two-earners probability distribu- 
tion function P2 (?') is 

P2(r)- Pi{r')Pi{r-r')dr' = ^e-^/''. (13) 

The function P2{r) il3i differs from the function Pi{r) by 
the prefactor r/R, which reflects the phase space available to 
compose a given total income out of two individual ones. It is 
shown as the solid curve in Fig.^] Unlike Pi (r), which has a 
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FIG. 15: Histogram: Probability distribution of income for families 
with two adults in 1996^^. Solid line: Fit to Eg. <13t . Inset histogram: 
Probability distribution of income for all families in 1996—. Inset 
solid line: 0.45Pi(r-) + 0.55P2(r). 



maximum at zero income, P2 (r) has a maximum at r = i? and 
looks qualitatively similar to the family income distribution 
curves in literature"*. 

From the same 1996 SIPP that we used in Sec. lIIlj^ \ we 
downloaded the variable TFTOTINC (the total family income 
for a month), which we then multiplied by 12 to get annual 
income. Using the number of family members (the variable 
EFNP) and the number of children under 18 (the variable 
RFNKIDS), we selected the families with two adults. Their 
distribution of family income is shown by the histogram in 
Fig. 01 The fit to the function il3\ . shown by the solid line, 
gives the parameter R listed in line (d) of TableU The families 
with two adults and more than two adults constitute 44% and 
1 1 % of all families in the studied set of data. The remaining 
45% are the families with one adult. Assuming that these two 
classes of families have two and one earners, we expect the 
income distribution for all families to be given by the super- 
position of Eqs. (|7} and O: 0.45Pi(r) + 0.55P2{r). It is 
shown by the solid line in the inset of Fig. ^] (with R from 
line (d) of TableU with the all families data histogram. 

By substituting Eq. il3\ into Eq. (jSji, we calculate the 
Lorenz curve for two-earners families: 

x{r) = 1 - (1 +f)e^^ (14) 
= x{f) (15) 

It is shown by the solid curve in Fig. [n] Given that x — 
y — f^exp(— f)/2 and dx = rcxp(— f)df, the Gini coef- 
ficient for two-earners families is: G2 = 2 {x — y) dx = 
exp{-2r) dr = 3/8 = 0.375. The points in Fig. [l7l 
show the Lorenz data and Gini coefficient for family income 
during 1947-1994 from Table 1 of Ref.^''. The Gini coeffi- 
cient is very close to the calculated value 0.375. 

The fundamental assumption which underlies the derivation 
of il3\ is the independence of the two incomes ri and r2 of 
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FIG. 16: Points: The income for families with two earners in 1999 
One family is represented by two points (r 1,7-2) and (r2,ri). The 
existence of correlations between ri and 7-2 is equivalent to a cluster- 
ing of points along the diagonal. The data shows little evidence for 
such correlation. 



the two-earner family. While the good fit of the function il3i 
is an implicit validation of this assumption, it would be good 
to show that the assumption holds true directly from data. The 
Census database does not provide the individual values of ri 
and r2, but only their sum r = ri + r2- We found that the 
PSID surve>*2i does have give for a family with two earners, 
each earners contribution to the total income of the family. We 
downloaded the variable ER16463 (total labor income, head) 
and ER16465 (total labor income, wife). We represent in Fig. 
[^each family by two points (ri,r2) and (r2,ri). The ex- 
istence of correlations between ri and r2 is equivalent to a 
clustering of points along the diagonal. Or stated otherwise, 
the independence of ri and r2 should result in a uniform dis- 
tribution of points on infinitesimal strips perpendicular to the 
diagonal. The plotted data shows little, if any evidence for 
such correlation. 

The distributions of the individual and family income dif- 
fer qualitatively. The former monotonically increases toward 
the low end and has a maximum at zero income (Fig. |6}. 
The latter, typically being a sum of two individual incomes, 
has a maximum at a finite income and vanishes at zero (Fig. 
I15> . Thus, the inequality of the family income distribution 
is smaller The Lorenz data for families follow the different 
Eq. jl5> . again without fitting parameters, and the Gini coef- 
ficient is close to the smaller calculated value 0.375 (Fig.llV^ 
Despite different definitions of income by different agencies, 
the parameters extracted from the fits (Table |lj are consistent, 
except for line (e). 

The qualitative difference between the individual and fam- 
ily income distributions was emphasized in Ref.^^, which split 
up joint tax returns of families into individual incomes and 
combined separately filed tax returns of married couples into 
family incomes. However, Refs.^^ and^^ counted only "indi- 
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FIG. 17: Solid curve: Lorenz plot <15> for distribution <13> . Points: 
Census data for families, 1947-1994"'. Inset points: Gini coefficient 
data for families from Census^'. Inset line: The calculated value 3/8 
of the Gini coefficient for distribution <13> . 



vidual tax returns", which also include joint tax returns. Since 
only a fraction of families file jointly, we assume that the latter 
contribution is small enough not to distort the tax returns dis- 
tribution from the individual income distribution significantly. 
Similarly, the definition of a family for the data shown in the 
inset of Fig. ^] includes single adults and one-adult families 
with children, which constitute 35% and 10% of all families. 
The former category is excluded from the definition of a fam- 
ily for the data^' shown in Fig.[nJ but the latter is included. 
Because the latter contribution is relatively small, we expect 
the family data in Fig.^]to approximately represent the two- 
earners distribution il3\ . Technically, even for the families 
with two (or more) adults shown in Fig.^J we do not know 
the exact number of earners. 

With all these complications, one should not expect perfect 
accuracy for our fits. There are deviations around zero income 
in Figs.|6lE] and^] The fits could be improved there by mul- 
tiplying the exponential function by a polynomial. However, 
the data may not be accurate at the low end because of under- 
reporting. For example, filing a tax return is not required for 
incomes below a certain threshold, which ranged in 1999 from 
$2,750 to $14,400^ As the Lorenz curves in Figs. Isl and [TTl 
show, there are also deviations at the high end, possibly where 
Pareto's power law is supposed to work. Nevertheless, the 
exponential law gives an overall good description of income 
distribution for the great majority of the population. 

Figs.|6land0demonstrate that the exponential law fits 
the individual income distribution very well. The Lorenz data 
for the individual income follow Eq. dlOl l without fitting pa- 
rameters, and the Gini coefficient is close to the calculated 
value 0.5 (Fig.©. 

It is interesting to study the distribution of the Gini index 
across the globe. We present in Fig. ^] data from a World 
Bank publication"*^. This is an question of utmost importance 
for the World Bank since the distribution of the Gini index 
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World distribution of Gini index, 1988 and 1993 
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FIG. 18: Distribution of Gini index for households across the globe 
for two different years, 1988 and 1992''** . For West Europe and North 
America, the Gini index is close to 0.375, the value predicted in Sec. 
Ill Fl A sharp increase in inequality is observed in the Eastern Europe 
and the Former Soviet Union (FSU) republics after the fall of the 
"iron curtain" when socialist economies were replaced by market- 
like economies. 



reflects tiie degree of inequality (poverty) in various regions 
of the globe. For West Europe and North America, the Gini 
index is close to 0.375, the value we predicted in Sec. Ill Fl A 
sharp increase in inequality is observed in the Eastern Europe 
and the Former Soviet Union (FSU) republics after the fall of 
the "iron curtain" when socialist economies where replaced 
by market-like economies. We conjecture that the high value 
of the Gini index in Asia may be due to the fact that a typical 
family in Asia has only one earner, so the Gini index is close 
to a value of 0.5 as observed for individuals. 

As the evidence from Fig.^jand Fig.[^shows the value of 
the Gini index varies in time and across the globe, but for the 
Western world it is close to the value of 0.375, as predicted in 
Sec.lHH 



G. Distribution of wealth 

In this section, we discuss the cumulative probability dis- 
tribution of wealth iV(ti;)=(the number of people whose in- 
dividual wealth is greater than w)/(the total number of peo- 
ple). A plot of N vs. w is equivalent to a plot of a person's 
rank in wealth vs. wealth, which is often used for top rich- 
est people-'^. We will use the power law, N{w) cx 1/w", 
and the exponential law N{w) cx exp{—w/W), to fit the 
data. These distributions are characterized by the exponent 
a and the "temperature" W. The corresponding probability 
densities, P{w) = —dN{w)/dw, also follow a power law or 
an exponential law. For the exponential law, it is also use- 
ful to define the temperatures VF(^) (also known as the me- 
dian) and 1^(10) using the bases of 1/2 and 1/10: N{w) cx 
(1/2)-/^*'' oc (1/10)-/^'"". 
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FIG. 19: Left panel: Cumulative probability distribution of to- 
tal net capital (wealth) shown in log-log, log-linear (inset) coordi- 
nates. Points: the actual data. Solid lines: fits to the exponential 
(Boltzmann-Gibbs) and power (Pareto) laws. 



The distribution of wealth is not easy to measure, because 
people do not report their total wealth routinely. However, 
when a person dies, all assets must be reported for the pur- 
pose of inheritance tax. Using these data and an adjustment 
procedure, the British tax agency, the Inland Revenue (IR), re- 
constructed wealth distribution of the whole UK population. 
In Figs. 119! and l20l we present the 1996 data obtained from 
their Web site^"^. Fig.ll9lshows the cumulative probability as 
a function of the personal total net capital (wealth), which is 
composed of assets (cash, stocks, property, household goods, 
etc.) and liabilities (mortgages and other debts). The main 
panel illustrates in the log-log scale that above 100 k£ the 
data follow a power law with the exponent a — 1.9. The 
inset shows in the log-linear scale that below 100 k£ the 
data is very well fitted by an exponential distribution with 
the temperature Wuk = 59.6 k£ (W^jl = 41.3 k£ and 



W, 



(10) 
UK 



= 137.2 ki;). 



Since we have estabilished that the distribution of wealth 
has an exponential regime followed by a power-law tail, we 
plot the Lorenz curve for wealth in the right panel of Fig. 1201 
As in Sec. lIIBI the horizontal and vertical coordinates are the 
cumulative population x{w) and the cumulative wealth y{w): 



x{w) 



P{w')dw\ 



(16) 



y{w)= / w'P[w')dw' / / w'P{w')dw'. (17) 
Jo Jo 

As in the case of US income data Fig.^| there is an abrupt 
change in the Lorenz curve for the top few percent of peo- 
ple. We interpret this fact in the same way as in Sec. IIIDI 
as a "Rose condensation of wealth". Again, we calculate / 
the fraction of the total wealth of the society if all people are 
described by the exponential law, to the actual total wealth 
of the society. For the United Kingdom in 1996, this frac- 



16 



United Kingdom, iR data for 1 996 United States, IRS data for 1 997 




FIG. 20: Points: Lorenz plot for wealth, United Kingdom 199&S. 
Solid curve: The Lorenz curve given by <18> with a condensate frac- 
tion of 6 = 16%. 



tion was / = 0.84, which means that the "condensate" has 
6=1 — / = 16% of the total wealth of the system. The 
functional form for the Lorenz curve for wealth is 

?/ = (1 - + (1 - x) ln(l - x)] + bS{l ~ x). (18) 

This function is plotted as a solid line in Fig. |20| One can 
see that the data systematically deviate from the exponential 
law because of the wealth concentrated in the power-law tail; 
however, the deviation is not very big. The so-called Gini 
coefficient'^'^^, which measures the inequality of wealth dis- 
tribution, has increased from 64% in 1984 to 68% in 1996^^ 
This value is bigger than the Gini coefficient 50% for a purely 
exponential distribution''^. 

While there seems to be no controversy about the fact that 
a few people hold a finite fraction of the entire wealth, the 
specific value of this fraction varies widely in the hterature. 
Fractions as high as 80-90% were reported in the literature^^ 
without any support from data. We provide a clear procedure 
of how to calculate this factor, and we find that b = 16%. 



H. Other distributions for income and wealtli 

As we discussed in the introduction. Sec. Ill Al there have 
been many proposed functional forms for the distribution of 
income and wealth. In this section we will investigate several 
of them. 

In Sees. Ill BllllDllllEI and HlGl we have shown ample evi- 
dence to support our findings that the distribution on income 
and wealth have a similar structure with an exponential part 
followed by a power-law tail. These results immediately in- 
validate the proposed Pareto or Gamma distributions as the 
correct distribution of income. 

Another popular distribution which has been proposed for 
many years as the distribution of income is the lognormal 



FIG. 21: Points: Calculated probability density P{r) for United 
States IRS income data for 1997 multiplied by income r, i.e. rP{r). 
Curve: Fit of data with function P{r) — rPLN{r) <19> with a 
quadratic function. Clearly, the lognormal distribution fails to de- 
scribe the data points accurately. 



distribution™. The lognormal distribution has the expression 



PLN{r) 



: exp 



(Inr — 



(19) 



where /i and a are two parameters. As given by (I19I I the 
functional form of the lognormal distribution makes it hard 
to compare it with data. But the function P{r) = rPiAr(r) 
is a quadratic function in a log-log plot. In Fig. |^ we plot- 
ted InP(r) versus Inr for the 1997 income data from IRS?^. 
As it can be seen from Fig.|2]a fit of P{r) with a quadratic 
function cannot be made. The data points show the power- 
law regime for high incomes. We conclude that the lognormal 
distribution does not describe the distribution of income. 

The distribution of wealth has essentially the same struc- 
ture as the distribution of income. An exponential regime 
for low wealth values and a power-law tail for high wealth 
values. Contrary to the distribution of income, the distribu- 
tion of wealth appears to have a continuous derivative in the 
cross-over between the two regimes (see the discussion at the 
beginning of Sec. Ill Dt . The shape of the wealth distribution 
suggests a modification of the regular exponential Boltzmann- 
Gibbs distribution for values of wealth w :s> W. In the past 
years, in the context of non-extensive statistical mechanics, 
Tsallis has proposed a viable alternative to the Boltzmann- 
Gibbs distribution'^*'. 

The Tsallis distribution can be obtained by maximizing 
a generalized entropy and is the subject of a lot of current 
research. The appeal of the Tsallis distribution is that it 
has power-law tails for large arguments, and that in certain 
limit it becomes the exponential distribution. The Tsallis 
distribution"**' has the expression 
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FIG. 22: Points: Internal Revenue wealth data for individuals, 
1996^^. Curve (a): Fit with the Kaniadakis distribution <21> . Curve 
(b): Fit with the Tsallis distribution <20> . Curve (c): Fit with the 
Bouchaud-Mezard/Solomon-Richmond distribution <22t . The inset 
shows the same three curves in a log-linear scale. The deviations of 
model <22t from data are evident for small values of income. 



where Wt is the dimensional parameter of the Tsallis distribu- 
tion, analogous to the temperature of the exponential distribu- 
tion. In the limit q 1, the distribution goes to the exponen- 
tial form Pq=i (w) cx q-'^'/'^t Pqj- g > the distribution has 
a power-law decay Pq{w) ^ {l/w)'^^'^~^ for values w ^ W. 

In Fig. |22] we present the fit of UK wealth data with the 
Tsallis distribution. Overall, the quality of the fit is good. The 
parameters of the Tsallis distribution implied by the fit are: 
q = 1.42, and for the Tsallis temperature Wt = 4:lk£. 

Another distribution proposed by Kaniadakis"*^ in the con- 
text of non-extensive statistical mechanics has the form of a 
deformed exponential 



1/k 



(21) 



where x = w/W^. For k ^ 0, the distribution has the ex- 
ponential form Pk=o{w) = e^""/^". For large arguments 
w ^ W, the distribution J21t has the power-law behavior 
Pniw) w (l/it;)^/". The fit with the Kaniadakis distribution 
is even better than that with the Tsallis distribution, as shown 
in Fig. |52 The parameters of the fit for the Kaniadakis dis- 
tribution are: k — 0.33 and for the Kaniadakis temperature 
= 45k£. 

A distribution for wealth put forward by Bouchaud and 
Mezard^^*, and independently by Solomon and Richmond^-, 
is 



Pc.{x) 



(a-l)"exp(-^) 
r(a) 



rct+l 



(22) 



where x = w/Wc- This distribution has an exponentially 
sharp cut-off for wealths w <C Wc, and a power-law tail 
Paiw) CX l/u;"+i for w > Wc- The fit of the UK wealth 



data with the distribution J22t is presented in Fig. |22| The 
values we obtain for the parameters are: a ~ 1.93 and 
Wc — 74 k£. The distribution J22t has its maximum at 
Xmax — {a ~ l)/(a + 1) = 0.32, which essentially states 
that for w > 0.32Wc — 23k£ the distribution is a power law. 
As seen in Fig.|22lthis claim is not supported by data. More- 
over, the percent of people in the power-law tail predicted by 
the distribution ( I22t . dxPa{x) — 81%, is clearly much 

larger than what is observed from the Lorenz curve for wealth, 

Fig.Hni ^^^^ 

All three proposed distributions ( 12012 11221 seem to capture 
well the power-law tail. In the inset of Fig.|22]a log-linear 
scale is used, and this makes the discrepancy between ( I22t 
and data much more evident. The distribution ( I22t is clearly 
inappropriate for low incomes. 



I. Conclusions 

Our analysis of the data shows that there are two clear 
regimes in the distribution of individual income. For low and 
moderate incomes up to approximately 95% of the total pop- 
ulation, the distribution is well described by an exponential, 
while the income of the top 5% individuals is described by a 
power-law (Pare to) regime. 

The exponential Boltzmann-Gibbs distribution naturally 
applies to quantities that obey a conservation law, such as en- 
ergy or money^^. However, there is no fundamental reason 
why the sum of incomes (unlike the sum of money) must be 
conserved. Indeed, income is a term in the time derivative 
of one's money balance (the other term is spending). Maybe 
incomes obey an approximate conservation law, or somehow 
the distribution of income is simply proportional to the distri- 
bution of money, which is exponential^^. 

Another explanation involves hierarchy. Groups of peo- 
ple have leaders, which have leaders of a higher order, and 
so on. The number of people decreases geometrically (ex- 
ponentially) with the hierarchical level. If individual income 
increases linearly with the hierarchical level, then the income 
distribution is exponential. However, if income increases mul- 
tiplicatively, then the distribution follows a power law^-. For 
moderate incomes below $100,000, the linear increase may 
be more realistic. A similar scenario is the Bernoulli trials^^, 
where individuals have a constant probability of increasing 
their income by a fixed amount. 

We found scaling in the cumulative probability distributions 
N{r) of individual income r derived from the tax statistics for 
different years in the UK and for different states in the US. The 
distributions Ni{r) have the scaling form f{r/Ri), where the 
scale Ri (the temperature) varies from one data set i to an- 
other, but the scaling function / does not. The function / has 
an exponential (Boltzmann-Gibbs) form at the low end, which 
covers about 95% of individuals. At the high end, it follows 
a power (Pareto) law with the exponents about 2. 1 for the UK 
and 1.7 for the US. Wealth distribution in the UK also has a 
qualitatively similar shape with the exponent about 1.9 and 
the temperature Wjjk = 60 k£. Some of the other values of 
the exponents found in literature are 1.5 proposed by Pareto 
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himself (a = 1.5), 1.36 found by Levy and Solomon for the 
distribution of wealth in the Forbes 400 list, and 2.05 found by 
Souma'*^ for the high end of income distribution in Japan. The 
latter study is similar to our work in the sense that it also uses 
tax statistics and explores the whole range of incomes, not just 
the high end. Souma^ finds that the probabiUty density P{r) 
at lower incomes follows a log-normal law with a maximum 
at a nonzero income. This is in contrast to our results, which 
suggest that the maximum of P{r) is at zero income. The dis- 
repancy may be due to the high threshold for tax reporting in 
Japan, which distorts the data at the low end. On the other 
hand, if the data is indeed valid, it may reflect the actual dif- 
ference between the social stuctures of the US/UK and Japan. 



III. DISTRIBUTION OF STOCK-PRICE FLUCTUATIONS 
A. Introduction 

Stochastic dynamics of stock prices is commonly described 
by a geometric (multiplicative) Brownian motion, which gives 
a log-normal probability distribution for stock price changes 
(returns)''". However, numerous observations show that the 
tails of the distribution decay slower than the log-normal 
distribution predicts (the so-called "fat-tails" effect)^ 
Particularly, much attention was devoted to the power-law 
tails^^-^^. The geometric Brownian motion model has two pa- 
rameters: the drift /i, which characterizes the average growth 
rate, and the volatility a, which characterizes the noisiness of 
the process. There is empirical evidence and a set of stylized 
facts indicating that volatility, instead of being a constant pa- 
rameter, is driven by a mean-reverting stochastic process^'^iS. 
Various mathematical models with stochastic volatility have 
been discussed in literatura^Si^Siffii^. 

In this chapter, we study a particular stochastic volatil- 
ity model, where the square root of the stock-price volatil- 
ity, called the variance, follows a random process known in 
financial literature as the Cox-lngersoll-Ross process and in 
mathematical statistics as the Feller process''^'^^. We solve the 
Fokker-Planck equation for this model exactly and find the 
joint probability distribution of returns and variance as a func- 
tion of time, conditional on the initial value of variance. The 
solution is obtained in two different ways: using the method of 
characteristics^^ and the method of path integrals AMj# . The 
latter is more familiar to physicists working in financ o^'i^^ . 

While returns are readily known from a financial time se- 
ries data, variance is not given directly, so it acts as a hidden 
stochastic variable. Thus, we integrate the joint probability 
distribution over variance and obtain the marginal probability 
distribution of returns wnconditional on variance. The latter 
distribution can be directly compared with financial data. We 
find excellent agreement between our results and the Dow- 
Jones data for the period of 1982-2001. Using only four fit- 
ting parameters, our equations very well reproduce the proba- 
bility distribution of returns for time lags between 1 and 250 
trading days. This is in contrast to popular ARCH, GARCH, 
EGARCH, TARCH, and similar models, where the number of 
parameters can easily go to a few dozenS. 



Our result for the probability distribution of returns has the 
form of a one-dimensional Fourier integral, which is easily 
calculated numerically or, in certain asymptotical limits, an- 
alytically. For large returns, we find that the probability dis- 
tribution is exponential in log-returns, which implies a power- 
law distribution for returns, and we calculate the time depen- 
dence of the corresponding exponents. In the limit of long 
times, the probability distribution exhibits scaling. It becomes 
a function of a single combination of return and time, with 
the scaling function expressed in terms of a Bessel function. 
The Dow-Jones data follow the predicted scaling function for 
seven orders of magnitude. 

The original theory of option pricing was developed 
by Black and Scholes for the geometric Brownian motion 
modeP". Numerous attempts to improve it using stochastic 
volatility models have been niade^^'^**-''^ ''"'*^ Particularly, op- 
tion pricing for the same model as in our pape^ was inves- 
tigated by Heston^". Empirical studies*^ show that Heston's 
theory fares better than the Black-Scholes model, but still 
does not fully capture the real-market option prices. Since 
our paper'*'^ gives a closed-form time-dependent expression 
for the probability distribution of returns that agrees with fi- 
nancial data, it can serve as a starting point for a better theory 
of option pricing. 

B. The model 

We consider a stock, whose price St, as a function of time t, 
obeys the stochastic differential equation of a geometric (mul- 
tiplicative) Brownian motion^: 

dSt = fiSt dt + atSt dW^f ^ • (23) 

Here the subscript t indicates time dependence, /i is the drift 
parameter, W^^^ is the standard random Wiener process^', and 
at is a time-dependent parameter, called the stock volatility, 
which characterizes the noisiness of the Wiener process. 

Since any solution of (I23> depends only on of, it is conve- 
nient to introduce the new variable vt = af, which is called 
the variance. We assume that vt obeys the following mean- 
reverting stochastic differential equation: 

dvt = ~-fivt-e)dt + Ky/^tdWt^^\ (24) 

Here 9 is the long-time mean of v, 7 is the rate of relaxation 
to this mean, Wt is the standard Wiener process, and k is a 
parameter that we call the variance noise. Eq. J24> is known 
in financial literature as the Cox-Ingersoll-Ross (CIR) process 
and in mathematical statistics as the Feller process''^ (p. 42). 
Alternative equations for vt, with the last term in (I24> replaced 
by K dWp^ or by nvt dWf'^\ are also discussed in literature^-. 
However, in this chapter, we study only the case given by Eq. 
(|24|i. 

We take the Wiener process appearing in i24i to be corre- 
lated with the Wiener process in ( I23t : 

dWt^^'' = pdWt^^^ + yjl- p^dZt, (25) 
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where Zt is a Wiener process independent of w\^^ , and p G 
[—1,1] is the correlation coefficient. A negative correlation 



(p < 0) between w'^'^ and Wl""' is known as the leverage 
effect" (p. 41). 

It is convenient to change the variable in ( I23t from price 
St to log-return rt — XvJyStj Sq). Using Ito's formula®, we 
obtain the equation satisfied by r* : 



.(2) 



dvt = (^p-j)dt + ^tdW^^\ 



(26) 



The parameter /i can be eliminated from i26\ by changing the 
variable to xt = rt — lit, which measures log-returns relative 
to the growth rate /i: 



dxt 



Vt 



dt + J¥tdW. 



(27) 



Where it does not cause confusion with r^, we use the term 
"log-return" also for the variable xt- 

Equations (I27> and \2A\ define a two-dimensional stochas- 
tic process for the variables xt and vt- This process is charac- 
terized by the transition probability Pt{x,v\vi) to have log- 
return X and variance v at time t given the initial log-return 
X = and variance Viait — 0. Time evolution of Pt (x, v\vi) 
is governed by the Fokker-Planck (or forward Kolmogorov) 
equation 



69 



(28) 



dx dv 



2dx^ 



The initial condition for ( I28> is a product of two delta func- 
tions 



Pt=o{x,v\ Vi) = 5{x) S{v - Vi). 



(29) 



The variance t; is a positive quantity, so Eq. (I28> is defined 
only for v > 0. However, when solving (I28> . it is convenient 
to extend the domain of v so that v e (— oo, oo). If 2'^9 > k^, 
such an extension does not change the solution of the equa- 
tion, because, given that P — for u < at the initial time 
t — 0, the condition Pt{x,v < 0) = is preserved for all 
later times t > 0. In order to demonstrate this, let us consider 
(l28i in the limit v — > 0: 



d_ 

dt' 



P 
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dP dP 

— +-lP + pn—. (30) 

ov ox 



Eq. ( I30> is a first-order partial differential equation (PDE), 
which describes propagation of P from negative v to positive 
V with the positive velocity ^9 — Thus, the nonzero 

function v > Q) does not propagate to w < 0, and 
Pt{x,v < 0) remains zero at all times t. Alternatively, it 
is possible to show^" (p. 67) that, if 2j9 > k^, the random 
process i24\ starting in the domain v > can never reach the 
domain t; < 0. 

The probability distribution of the variance itself, 
Ilt{v) — J dx Pt{x,v), satisfies the equation 



d_ 
di 



nt{v) = ^ Mv - 9)Ut{v)] + , (31) 



which is obtained from ( I28t by integration over x. Eq. (13 U 
has the stationary solution 



n,(v) 
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(32) 

which is the Gamma distribution. The maximum of n*(ti) is 
reached at v^ax — (3/a — 9 ~ The width w of n*(ti) 

can be estimated using the curvature at the maximum w ~ 
{k^ 1 2-^1 ^2-^9 1 - 1. The shape of YL^(y) is characterized 
by the dimensionless ratio 



276I 



(33) 



When 2^9/ 00, x ^ 00 and n*(?;) 5{v 



C. Solution of the Fokker-Planck equation 



Since x appears in ( I28t only in the derivative operator 
d/dx, it is convenient to make the Fourier transform 

1 r+oo 

P,{x,v\v,)^— dp,^P^^Pt,pAv\v^). (34) 
Inserting (I34t into \2%\ . we find 

lf=^i[i^-m (35) 

V — IPKP^-—V TT^V 

2 ^ ^""dv 2 dv^ 



Eq. i35\ is simpler than ( I28> . because the number of variables 
has been reduced to two, v and t, whereas only plays the 
role of a parameter. 

Since Eq. i35\ is linear in v and quadratic in d/dv, it can 
be simplified by taking the Fourier transform over v 

Pt,p. iv\v,) = ^ r dp, e'P^^Pt,p^ ip, \v,). (36) 

The PDE satisfied by Pt,p^ {pv \ Vi) is of the first order 

d 



d-t + V^^^-p^ + ^^j dp. 



P = -i-t9p.,P, 
(37) 



where we introduced the notation 

T = "f + ipKpx. (38) 
Eq. ( I37> has to be solved with the initial condition 

Pt=o,p^ ipv \vt) = exp(-ipt,t)j). (39) 

The solution of the PDE (j^J is given by the method of 
characteristics''^: 



Pt.p^ [Pv I Vi) = exp ( -ipy{0)v^ - ij9 / drp^ir) 

Jo J 



(40) 
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where the function p„ (r) is the solution of the characteristic 
(ordinary) differential equation 

J^^rp4T) + —pl{T) + -{pl-2p,) (41) 

with the boundary condition p„ (t) = p^, specified at t t. 
The solution ( I40> can be also obtained using the method of 
path integrals described in Sec. 1111 Dl The differential equa- 
tion ( I41> is of the Riccati type with constant coefficients^^, 
and its solution is 



Pv{t) = -i 



2n 1 



+ i- 



.T-n 



where we introduced the frequency 



n = ^T^ + K^ipl-ip,). 

and the coefficient 

2n 



K^p„ - i{V - Vl) ' 



(42) 



(43) 



(44) 



Substituting J42> into J40> and taking the Fourier trans- 
forms (I34> and i36\ . we get the solution 



Pt{x,v\vi) 



X exp < —ipy{Q)vi + 



dpx dpy e 



je{T~n)t 276I c-e 



(45) 



In- 



of the original Fokker-Planck equation ( I28l l with the initial 
condition ( I29> . where Pt,(T = 0) is given by ( I42> . 



D. Path-Integral Solution 

The Fokker-Planck equation ( I35> can be though of as a 
Schrodinger equation in imaginary (Euclidean) time: 

d— . _ 

■Q-^Pt,pAv\vi) = ~Hp^{p^,v)Pt,p^{y\vi) (46) 

with the Hamiltonian 

^ ^ p %p 

H = -^Pv^ ~ ^IP-^i^ ^0) + — + pupxPvV. (47) 

In ( I47> we treat pu and ?) as canonically conjugated op- 
erators with the commutation relation [v,Pv] — i- The 
transition probability Pp_^{v,t\vi) is the matrix element of 
the evolution operator exp(— iJt) and has a path-integral 
representationSi^i^ 

(48) 

Here the action functional Sp^ [Pv{t), v{t)] is 

Sp^ = / dT {tp,{T)i{T) - Hp^ [p,{t),v{t)]} , (49) 
Jo 



and the dot denotes the time derivative. The phase-space path 
integral ( I48t is the sum over all possible paths Pv{t) and v{t) 
with the boundary conditions v{t = Q) ~ Vi and v{t — t) = 
V imposed on v. 

It is convenient to integrate the first term on the r.h.s. of ( I49> 
by parts: 



Sp:, = i[Pv{t)v ~ Py{0)Vi 



- / dT 

Jo 



5H 



Sv(t) 



- 176* / dTPy{T) 

JO 

v{t), 



(50) 



where we also separated the terms linear in v{t) from the 
Hamiltonian (|^}. Because v(r) enters linearly in the action 
(15 0> . taking the path integral over Vv generates the delta- 
functional 5 [pv{t) — Pv{t)], where p„(t) is the solution of 
the ordinary differential equation ( I41> with a boundary con- 
dition specified at r = i. Taking the path integral over Vp^ 
resolves the delta-functional, and we find 



/ I \ f dpy y—p (Q\y-'\—-i'^g C drp^fr) 

,pj.v\vt) = I ^'^ ^ ' ' ^ 'Jo '\ 



(51) 



where p„(t = t) = Py. Eq. <51l l coincides with j40> after the 
Fourier transform J36l l. 



E. Averaging over variance 

Normally we are interested only in log-returns x and do not 
care about variance v. Moreover, whereas log-returns are di- 
rectly known from financial data, variance is a hidden stochas- 
tic variable that has to be estimated. Inevitably, such an esti- 
mation is done with some degree of uncertainty, which pre- 
cludes a clear-cut direct comparison between v\vi) and 
financial data. Thus, we introduce the probability distribution 



+00 



Pt{x \vi)= / dvPt{x,v\ V,), 



(52) 



where the hidden variable v is integrated out. The integration 
of ( I45> over v generates the delta-function 6{py), which ef- 
fectively sets py = 0. Substituting the coefficient C from i44\ 
with Py = into (I45> . we find 



Ptix\v^) 



1 

2^ 



+ 00 



dpaje'^""" ''■r+nooth(nt/2) 



X e 



■ ln(cosh ^+s7 sinh 



(53) 



To check the validity of ( I53> . let us consider the case k = 0. 
In this case, the stochastic term in (I24t is absent, so the time 
evolution of variance is deterministic: 

v{t) = + {v, - 0)e-''K (54) 

Then process ( I27> gives a Gaussian distribution for x, 



>(k=0) 



(X I Vi) = 



: exp 



{x + vtt/2f 
2vtt 



(55) 
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with the time-averaged variance Vt = j Jq dTv{T). Eq. (I55> 
demonstrates that the probability distribution of stock price S 
is log-normal in the case k = 0. On the other hand, by taking 
the limit k ^ Q and integrating over in ( I53t . we reproduce 
the same expression i55\ . 

Eq. i53\ cannot be directly compared with financial time 
series data, because it depends on the unknown initial variance 
Vi. In order to resolve this problem, we assume that Vi has the 
stationary probability distribution 11* (ui), which is given by 
j32L^^ Thus we introduce the probability distribution Pt{x) 
by averaging ( I53I I over Vi with the weight 11, (v^): 



Dow Jones data, 1982-2001 



Pt{x) = / dv,Il4v.,)Pt{x\v^). 



(56) 



The integral over Vi is similar to the one of the Gamma func- 
tion and can be taken explicitly. The final result is the Fourier 
integral 



1 

2^ 



-|-oo 



(57) 



with 



Ft {p.) 
276* 



In 



cosh ■ 



nt 02 _ p2 



2717 



27r . , nt 

Sinn — 

2 



(58) 



The variable enters ( I58I I via the variables F from ( I38I I and 
from ( I43> . It is easy to check that Pt{x) is real, because 
ReF is an even function of p^ and Imi^ is an odd one. One 
can also check that Ft {px = 0) = 0, which implies that Pt (x) 
is correctly normalized at all times: J dx Pt{x) — 1. The 
second term in the r.h.s. of ( I58> vanishes when p — 0, i.e. 
when there are no correlations between stock price and vari- 
ance. The simplified result for the case p — Q is given in llll Fl 
by Eqs. O, and (El. 

Eqs. (I57> and (I58> for the probability distribution Pt{x) of 
log-return x at time t are the central analytical result of the 
chapter. The integral in (|^} can be calculated numerically 
or, in certain regimes discussed in Sees. 1111 Fl and UlIGl an- 
alytically. In Fig. |23l the calculated function Pt{x), shown 
by solid lines, is compared with the Dow- Jones data, shown 
by dots. (Technical details of the data analysis are discussed 
in Sec. 1111 Ffl ) Fig. 1231 demonstrates that, with a fixed set of 
the parameters 7, 9, k, p., and p, Eqs. i57i and (I58> very well 
reproduce the distribution of log-returns x of the Dow-Jones 
index for all times t. In the log-linear scale of Fig.|23] the tails 
of In Pt{x) vs. X are straight lines, which means that that tails 
of Pt (x) are exponential in x. For short times, the distribution 
is narrow, and the slopes of the tails are nearly vertical. As 
time progresses, the distribution broadens and flattens. 



F. Asymptotic behavior for long time t 

Eq. i24l implies that variance reverts to the equilibrium 
value 9 within the characteristic relaxation time I/7. In this 
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FIG. 23: Probability distribution Pt{x) of log-return x for different 
time lags t. Points: The Dow-Jones data for t — 1,5, 20, 40, and 250 
trading days. Solid lines: Fit of the data with Eqs. <57t and j58> . For 
clarity, the data points and the curves for successive t are shifted up 
by the factor of 10 each. The inset shows the curves without vertical 
shift. 



section, we consider the asymptotic limit where time t is much 
longer than the relaxation time: 7< ^ 2. According to ( 13 8> 
and ( I43> . this condition also implies that fit ^ 2. Then Eq. 
dSU reduces to 



Ft{p.)^^{T-n). 

Let us change of the variable of integration in ( I57> to 



Px 



■ Px +ipa, 



(59) 



(60) 



where 



Substituting (jSO) into (j^H), (ED, and (|55), we transform ( BTt 
to the following form 



where 

A = 

and 



-pox+At roo 



dp, cos(Ap,)e-^Vi+^, (62) 



j9t 



x + p— , B 



j9ujQt 



"f9 27 — pK 

2k2 l-p2 ■ 



(63) 



(64) 



According toS, the integral in (|62} is equal to 
BKiiVA^ + B^)/\/A^ + 52, where Ki is the first-order 
modified Bessel function. 
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Dow-Jones data, 1982-2001 
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FIG. 24: Renormalized probability density Pt{x)e'^°'^ / Nt plotted as 
a function of the scaling argument z given by Solid line: The 
scaling function P* (z) = isTi (2) from <65> . where Ki is the first- 
order modified Bessel function. Symbols: The Dow-Jones data for 
different time lags t. 
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FIG. 25: The fraction Gt of the total probability contained in the 
Gaussian part of Pt{x) vs. time lag t. Inset: Time dependence of the 
probability density at maximum Pt{xm) (points), compared with the 
Gaussian t"^/^ behavior (solid line). 



Thus, Eq. \51\ in the hmit 7^ ^ 2 can be represented in the 
scahng form 

Pt{x) = Nte-P°^P4z), P^z) = Ki{z)/z, (65) 
where the argument z = \/ + B"^ is 



{x + fr/9t/Ky 

1 



and the time-dependent normalization factor Nt is 



(66) 



(67) 



Eq. j65> demonstrates that, up to the factors Nt and e"^'''^, the 
dependence of Pt {x) on the two arguments x and t is given by 
the function P*(z) of the single scaling argument z in \66\ . 
Thus, when plotted as a function of z, the data for different 
X and t should collapse on the single universal curve P, (z). 
This is beautifully illustrated by Fig.|23 where the Dow- Jones 
data for different time lags t follows the curve (z) for seven 
decades. 

In the limit z ^ 1, we can use the asymptotic expression^^ 
Ki{z) « e~^A/7r/2z in (l65i and take the logarithm of P. 
Keeping only the leading term proportional to z and omitting 
the subleading term proportional to In z, we find 



In 



Nt 



-pox — z for z 3> 1. 



(68) 



Let us examine Eq. ( I68> for large and small 

In the first case \x\ ^ jOt/n, Eq. (I66> gives z 
u)q\x\I K^f\^^(? , SO Eq. ( I68> becomes 



In 



Nt 



-PqX 



^0 



\x\. 



(69) 



Thus, the probability distribution Pt{x) has the exponential 
tails (I69> for large log-returns \x\. Notice that, in the consid- 
ered limit 7^ ^ 2, the slopes d\nP/dx of the exponential 
tails (I69> do not depend on time t. Because of pq, the slopes 
\69\ for positive and negative x are not equal, thus the distri- 
bution Pt{x) is not symmetric with respect to up and down 
price changes. According to ( I61> , this asymmetry is enhanced 
by a negative correlation p < between stock price and vari- 
ance. 

In the second case l^j ^ ^6t/ n, by Taylor-expanding z in 
(I66> near its minimum and substituting the result into (16 8> . we 
get 



In 



Pt{x) 



-Pox 



2(1- p^)j0t 



(70) 



where Nj. — Ntexp{—u!o'^9t/ k'^). Thus, for small log- 
returns |a;|, the probability distribution Pt{x) is Gaussian with 
the width increasing linearly in time. The maximum of Pt{x) 
in ( 17 0> is achieved at 



p{luo - 7) 



(71) 



Eq. (I71> gives the most probable log-return Xm (t) at time t, 
and the coefficient in front of t constitutes a correction to the 
average growth rate so that the actual growth rate is p. = 

fi + dxm/dt. 

As Fig.|23lillustrates, In Pt (x) is indeed linear in x for large 
\x\ and quadratic for small in agreement with ( I69t and 
( 17 Ot . As time progresses, the distribution, which has the scal- 
ing form i65i and ( I66t . broadens. Thus, the fraction Gt of the 
total probability contained in the parabolic (Gaussian) portion 
of the curve increases, as illustrated in Fig. 1251 

In the remainder of the section, we explain the way the 
Gaussian weight of the distribution has been calculated. As 
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explained in Sec. 1111 HI the case relevant for comparison with 
the data is p = 0. In this case, by shifting the contour of 
integration in i57\ as follows px ^ Px + we find 



where 

FtiPx 
and 



72 61^ 276I 



In 



2tt 



cosh 1 smn — 

2 27f7 2 



n=^-i^ + k^pI + 1/4). 



(73) 



(74) 



Now the function Ft{px) is real and symmetric: Ft{px) — 
Ft {—Px)- Thus, the integral in M2\ is a symmetric function of 
X. So it is clear that the only source of asymmetry of Pt {x) in 
X is the exponential prefactor in Mil , as discussed at the end 
of Sec. IIII Hi Eqs. (I72> . (I73> . and MA\ are much simpler than 
those for p 7^ 0. 

Let us expand the integral in M2\ for small x: 



(75) 



where the coefficients are the first and the second moments of 

exp[Ft(p^)] 



+ 00 



dpx 
27r 



2tt 



(76) 

On the other hand, we know that Pt (x) is Gaussian for small 
X. So we can write 



(77) 



with the same coefficients as in ( I75> . If we ignore the exis- 
tence of fat tails and extrapolate to x E (—00, 00), the 
total probability contained in such a Gaussian extrapolation 
will be 



Gt= / da;/ioe-"/2"^^^'/2M« ^ W^e'^«/«'^^ (78) 



Obviously, Gt < 1, because the integral ( I78> does not take 
into account the probability contained in the fat tails. Thus, 
the difference 1 — Gt can be taken as a measure of how much 
the actual distribution Pt{x) deviates from a Gaussian func- 
tion. 

We calculate the moments i76\ numerically for the function 
F given by ( I73> . then determine the Gaussian weight Gt from 
( I78> and plot it in Fig.|25]as a function of time. For t 00, 
Gt ~* 1, i.e. Pt{x) becomes Gaussian for very long time lags, 
which is known in literature^^. In the opposite limit < ^ 0, 
Ft{px) becomes a very broad function of px, so we cannot 
calculate the moments po and /12 numerically. The singular 
limit t ^ requires an analytical study. 
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FIG. 26: Complex plane of px- Dots: The singularities of Ftipx)- 
Circled crosses: The accumulation points ±iq^ of the singularities 
in the limit 7* ^ 2. Symbol x : Saddle point ps, which is located 
in the upper half-plane for x > Q. Dashed line: The contour of 
integration displaced from the real axis in order to pass through the 
saddle point ps. 



Fig.|25]shows that, at sufficiently long times, the total prob- 
ability contained in the non-Gaussian tails becomes negligi- 
ble, which is known in literature^^. The inset in Fig.l25lillus- 
trates that the time dependence of the probability density at 
maximum, Pt{xm)^ is close to t^^/^, which is characteristic 
of a Gaussian evolution. 



G. Asymptotic behavior for large log-return x 

In the complex plane of px, function F{px) becomes sin- 
gular at the points px where the argument of any logarithm 
in ( I58> vanishes. These points are located on the imaginary 
axis of Px and are shown by dots in Fig. The singu- 
larity closest to the real axis is located on the positive (neg- 
ative) imaginary axis at the point p'l {pi). At these two 
points, the argument of the last logarithm in ( I58t vanishes, and 
we can approximate F{px) by the dominant, singular term: 
F{px)^-{2-ie/K^)\n{px-p^). 

For large |a:|, the integrand of \51\ oscillates very fast as a 
function of px- Thus, we can evaluate the integral using the 
method of stationary phase^*^ by shifting the contour of inte- 
gration so that is passes through a saddle point of the argument 
ipxX + F{px) of the exponent in i5H . The saddle point posi- 
tion Ps, shown in Fig.|26|by the symbol x, is determined by 
the equation 



dF{px) 



dpx 



276I 
T2" 



P^=Ps 



x>Q, 
P.-P. (79) 

, a; < 0. 



Ps-P 

For a large |a;| such that \xpf \ ^ 276'/^^, the saddle point 
is very close to the singularity point: ps ~ p^ for x > Q and 
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FIG. 27: small Solid line: The slope = — d In P/ dx of the expo- 
nential tail for a:: > as a function of time. Points: The asymptotic 
approximation <82> for the slope in the limit <^ 2. Dashed line: 
The saturation value qt for 7t 3> 2, Eq. <81> . 
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FIG. 28: Log-linear plot of the Dow- Jones versus time, for the pe- 
riod 1982-2001. The straight solid line through the data points repre- 
sents a fit of the Dow- Jones index with an exponential function with 
growth rate p = 13.3%/year. 



Ps ~ Pi for a; < 0. Then the asymptotic expression for the 
probability distribution is 



H. Comparison with Dow- Jones time series 



, x> 0, 
x < 0, 



(80) 



where qf = T'ipfi't) ^"e real and positive. Eq. ( I80> shows 
that, for all times t, the tails of the probability distribution 
Pt{x) for large |a;| are exponential. The slopes of the exponen- 
tial tails, — =F d In P/dx, are determined by the positions 
pf of the singularities closest to the real axis. 

These positions pf{t) and, thus, the slopes depend 
on time t. For times much shorter than the relaxation time 
(7t <C 2), the singularities lie far away from the real axis. 
As time increases, the singularities move along the imaginary 
axis toward the real axis. Finally, for times much longer than 
the relaxation time {'ft ^ 2), the singularities approach lim- 
iting points: pf ^iqf^ which are shown in Fig. 1261 bv 
circled x 's. Thus, as illustrated in Fig. ^\ the slopes qf 
monotonously decrease in time and saturate at long times: 



Qt 



±Pa 



LUq 



for 7i > 2. 



(81) 



The slopes J81> are in agreement with Eq. (I69> valid for 7^ ^ 
2. The time dependence qf at short times can be also found 
analytically: 
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K2(l-p2)i 



for 7t < 2. (82) 



The dotted curve in Fig. ^\ shows that Eq. ( I82> works very 
well for short times t, where the slope diverges at i ^ 0. 



To test the model against financial data, we downloaded 
daily closing values of the Dow-Jones industrial index for the 
period of 20 years from 1 January 1982 to 31 December 2001 
from the Web site of Yahoc^^. The data set contains 5049 
points, which form the time series {Sr}, where the integer 
time variable r is the trading day number. We do not filter the 
data for short days, such as those before holidays. 

Given {St}, we use the following procedure to extract the 
probability density p^^'^^ {r) of log-return r for a given time 
lag t. For the fixed t, we calculate the set of log-returns 
— \nSr+t/ Sr} for all possible times r. Then we par- 
tition the r-axis into equally spaced bins of the width Ar and 
count the number of log-returns belonging to each bin. In 
this process, we omit the bins with occupation numbers less 
than five, because we consider such a small statistics unreli- 
able. Only less than 1 % of the entire data set is omitted in this 
procedure. Dividing the occupation number of each bin by 
Ar and by the total occupation number of all bins, we obtain 



the probability density P^ 
find pj:^'^\x), we replace 



(r) for a given time lag t. To 
+ /it. In Fig.|28lwe plot 



the value of the Dow- Jones index, and the fit with an expo- 
nential function corresponding to an inflation rate ^. In Fig. 
1291 we illustrate the histogram corresponding to the probabil- 
ity density p^^''^ (r) for a time lag t = 5 days. The first and 
last bin in Fig. |29] show all the events outside the two vertical 
dashed lines. These events are very rare, contribute less than 
1% to the entire number of data points, and therefore they are 
discarded in the fitting process. 

Assuming that the system is ergodic, so that ensemble aver- 
aging is equivalent to time averaging, we compare p^*-^"^-* (x) 
extracted from the time series data and Pt{x) calculated in 
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Dow-Jones data, 1982-2001 

500 1 ^ ^ ^ : 




FIG. 29: Histogram of five day log-retums rt^s. The first and the 
last bin, which are separated from the rest by the two vertical dashed 
lines, are discarded in the fitting process. The population of the 
first (last) bin contains all the extreme events from — co(+cxd) to the 
left(right) dashed line. 



previous sections, which describes ensemble distribution. In 
the language of mathematical statistics, we compare our the- 
oretically derived population distribution with the sample dis- 
tribution extracted from the time series data. We determine 
parameters of the model by minimizing the mean-square devi- 
ation J In P/^-^^ (x) - In Pt(a;)|2, where the sum is taken 
over all available x and t = 1, 5, 20, 40, and 250 days. 
These values of t are selected because they represent differ- 
ent regimes: jt 1 for t ~ 1 and 5 days, jt ^ 1 for t — 20 
days, and •jt ^ 1 for t = 40 and 250 days. As Figs.l23land 
illustrate, our expression ( I57l i and ( I58l l for the probability 
density Pt{x) agrees with the data very well, not only for the 
selected five values of time t, but for the whole time interval 
from 1 to 250 trading days. However, we cannot extend this 
comparison to t longer than 250 days, which is approximately 
1/20 of the entire range of the data set, because we cannot 
reliably extract p^^'^^ (^x) from the data when t is too long. 

The values obtained for the four fitting parameters (7, 6, k, 
/J,) are given in Tablelllll Within the scattering of the data, we 
do not find any discernible difference between the fits with the 
correlation coefficient p being zero or slightly different from 
zero. Thus, we conclude that the correlation p between the 
noise terms for stock price and variance in Eq. ( I25t is practi- 
cally zero, if it exists at all. Our conclusion is in contrast with 
the value p = —0.58 found in^"* by fitting the leverage coiTe- 
lation function introduced in^^. Further study is necessary in 
order to resolve this discrepancy. Nevertheless, all theoretical 
curves shown in this chapter are calculated for p — 0, and they 
fit the data very well. 

All four parameters (7, 9, k, p) shown in Tablelllllhave the 
dimensionality of 1/time. The first line of the Table gives their 
values in the units of 1/day, as originally determined in our fit. 
The second line shows the annualized values of the parameters 



Units 


7 


e 


K 




I/day 


4.50 X 10"^ 


8.62 X 10"^ 


2.45 X 10"^ 


5.67 X 10"'' 


1/year 


11.35 


0.022 


0.618 


0.143 



TABLE III: Parameters of the model obtained from the fit of the 
Dow- Jones data. We also find p = for the correlation coefficient 
and 1/7 — 22.2 trading days for the relaxation time of variance. 



in the units of 1/year, where we utilize the average number of 
252.5 trading days per calendar year to make the conversion. 
The relaxation time of variance is equal to I/7 = 22.2 trading 
days = 4.4 weeks « 1 month, where we took into account that 
1 week = 5 trading days. Thus, we find that variance has a 
rather long relaxation time, of the order of one month, which 
is in agreement with the conclusion of Ref.^"*. 

Using the numbers given in Table |in| we find the value 
of the parameter 276'/k^ ~ 1.296. Since this parameter is 
greater than one, the stochastic process ( I24> never reaches the 
domain t; < 0, as discussed in Sec. lIIIBl 

The effective growth rate of stock prices is determined by 
the coordinate r„i{t) where the probability density Pt{rm) is 
maximal. Using the relation — Xm + pt and Eq. ( 17 1> . we 
find that the actual growth rate is p, = p — 76/2^0 — 13% per 
year. This number coincides with the average growth rate of 
the Dow -Jones index obtained by a simple fit of the time series 
{Sr} with an exponential function of r. The effective stock 
growth rate p is comparable with the average stock volatility 
after one year a = \fB — 14.7%. Moreover, the parameter 
( I33> . which characterizes the width of the stationary distribu- 
tion of variance, is equal to x = 0.54. This means that the 
distribution of variance is broad, and variance can easily fluc- 
tuate to a value twice greater than the average value Q. Thus, 
even though the average growth rate of stock index is posi- 
tive, there is a substantial probability dr Pt(r) ~ 17.7% 
to have negative growth for t = 1 year. 

According to j81> . the asymmetry between the slopes of 
exponential tails for positive and negative x is given by the 
parameter po> which is equal to 1/2 when p = (see also the 
discussion of Eq. ( I72> at the end of Sec. IIII F> . The origin 
of this asymmetry can be traced back to the transformation 
from (l23t to (l26t using Ito's formula. It produces the term 
0.5uf dt in the r.h.s. of (I26> . which then generates the second 
term in the r.h.s. of ( I28> . The latter term is the only source of 
asymmetry in x of Pt {x) when p — 0. However, in practice, 
the asymmetry of the slopes po = 1/2 is quite small (about 
2.7%) comparted with the average slope loq/h— 18.4. 



I. Conclusions 

We derived an analytical solution for the probability dis- 
tribution Pt{x) of log-returns x as a function of time t for 
the model of a geometrical Brownian motion with stochastic 
variance. The final result has the form of a one-dimensional 
Fourier integral (|^} and ( I58> . [In the case p ~ Q, the 
equations have the simpler form M2\ , ( I73> . and (I74M Nu- 
merical evaluation of the integral (I57> is simple compared 
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with computationally-intensive numerical solution of the orig- 
inal Fokker-Planck PDE or Monte-Carlo simulation of the 
stochastic processes. 

Our result agrees very well with the Dow-Jones data, as 
shown in Fig.|23] Comparing the theory and the data, we de- 
termine the four (non-zero) fitting parameters of the model, 
particularly the variance relaxation time I/7 = 22.2 days. 
For time longer than I/7, our theory predicts scaling behavior 
i65\ and i66\ . which the Dow-Jones data indeed exhibits over 
seven orders of magnitude, as shown in Fig. |23 The scal- 
ing function P^,{z) = Ki{z)/z is expressed in terms of the 
first-order modified Bessel function Ki. Previous estimates 
in literature of the relaxation time of volatility using various 
indirect indicators range from 1.5 days^^ (p. 80) to less than 
one day and more than few tens of days^' (p. 70) for S&P 500 
and 73 days for the half-life of the Dow-Jones index^''. Since 
we have a very good fit of the entire probability distribution 
function for times from 1 to 250 trading days, we believe that 
our estimate, 22.2 days, is much more reliable. A close value 
of 19.6 days was found in Refill. 

As Fig. 1231 shows, the probability distribution Pt{x) is ex- 
ponential in X for large |a;|, where it is characterized by time- 
dependent slopes d\nP/dx. The theoretical analysis pre- 
sented in Sec. IIIICil shows that the slopes are determined by 
the singularities of the function Ft{px) from ( I58> in the com- 
plex plane of that are closest to the real axis. The cal- 
culated time dependence of the slopes d\nP/dx, shown in 
Fig.|^ agrees with the data very well, which further supports 
our statement that 1/7 = 22.2 days. Exponential tails in the 
probability distribution of stock log-returns have been noticed 
in literature befora^ (p. 61)p^ however time dependence of 
the slopes has not been recognized and analyzed theoretically. 
In Ref.^^, the power-law behavior of the tails has been em- 
phasized. However, the data for S&P 500 were analyzed in 
Ref.^^ only for short time lags t, typically shorter than one 
day. On the other hand, our data analysis is performed for the 
time lags longer than one day, so the results cannot be directly 
compared. 

Our analytical expression for the probability distribution 
of returns can be utilized to calculate option pricing. No- 
tice that it is not necessary to introduce the ad-hoc function 
X{S,v,t), the market-price of risk, as is typically done in 
literatura^SiSSSiSiffi, because now Pt (x) is known explicitly. 
Using the true probability distribution of returns Pt (x) should 
improve option pricing compared with the previous efforts in 
literature. 

Although we tested our model for the Dow-Jones index, 
there is nothing specific in the model which indicates that it 
applies only to stock market data. It would be interesting to 
see how the model performs when applied to other time series, 
for example, the foreign exchange dataS, which also seem to 
exhibit exponential tails. 



IV. PATH-INTEGRAL SOLUTION OF THE 
COX-INGERSOLL-ROSS/FELLER MODEL 

In this section we provide details on the path-integral solu- 
tion of the variance process (I24> . The stochastic differential 
equation of the Cox-Ingersoll-Ross/Feller process, Eq. J24> . 
is 

dxt = -j{xt - e)dt + K^tdWt. (83) 

Here 6 is the long-time mean of x, 7 is the rate of relaxation 
to this mean, Wt is the standard Wiener process, and k is a 
parameter that controls the amplitude of the noise. Eq. (18 3> is 
known in financial literature as the Cox-Ingersoll-Ross (CIR) 
process and it was used to model of interest rate dynamics. 
Feller- 2 used equation (I83> to model population extinction in 
biological systems. Recently, the same equation ( I83> was used 
to model the thermal reversal of magnetization in a magnetic 

• 80 

gram . 

One quantity of interest to calculate is the conditional (tran- 
sition) probability P{x, t \ Xi, ti), which gives the probability 
to find the system at position x at time t, given that it was 
at position Xi at initial time ti. The transition probability 
P{x, 1 1 Xi,ti) of the stochastic process described by (I83> sat- 
isfies the Fokker-Planck equation 

d d 

—Pix, 1 1 x„ U) = — [7(x - 6)P{x, 1 1 x,,ti)] 

^2 q2 

+ —^[xP{x,t\x^,U)] (84) 
2 o'^x 

The Fokker-Planck equation ( I84> can be though of as a 
Schrodinger equation in imaginary (Euclidean) time written 
in the coordinate "x" representation, 

t\x„U)^ ~H{p, x)P{x, 1 1 x„U) (85) 

ot 

with the Hamiltonian 

H = —fx-i-ip{x-e). (86) 

In the coordinate representation, the momentum operator p is 
given by the derivative operator p — —id/dx. The canoni- 
cal commutation relation holds true: [x,p] = i. Using this 
analogy with quantum mechanics, the transition probability 
P{x, 1 1 Xi,ti) can be found using the method of path inte- 
grals. The transition probability is the matrix element of the 
evolution operator U{t,ti) ~ exp{—H{t — ti)) between the 
initial state \xi), and final state {xf \, 

Pixf,tf\x,,t,) = {xf\e-"^'f-'^'>\x,). (87) 

The standard construction of the path integral is obtained 
by splitting the time interval [ti,tf] into M time steps of equal 
size e = {tf — ti)/M. One uses the multiplication property of 
the transition probability to write 

» Af-l 
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M 



(88) 



The contribution coming from interval {tj, ij+i) is evaluated 

as 



-H(t,-t,_i) 



(89) 



The matrix element inside the integral in Eq.( I89> can be now 
evaluated to 

(p,|e-^(*^-*-i)|x,_i) « e-^^(f-^-i)(p,>j_i). (90) 
Substitute the scalar product 



{Xj\pj 



'2tt 



{Vj\xj-i) = 



2tt 



(91) 



and Eq.(|90j back into (|89}, to obtain 



P{Xf,tf\Xi,ti) 



M-l . M , 



where 



H{pj,Xj^i) = —p^Xj^i - hPj{xj_i 



(92) 



(93) 



By taking the limit M 
integral expression 



(94) 



oo in Eq.(l92li. one arrives at the path 

Pixf,tf\x,,U) = JvxVpe'^^^'P^ 
where the action functional S[xt,pt] is given by 

S[xt,pt] 



dt 



iptxt - -l^Pt^t + ilPt [xt - 0) 



(95) 

The phase-space path integral from Eq.(l94li has to be calcu- 
lated over all continuous paths p{t) and x{t) satisfying the 
boundary conditions x{t,i) = Xi and x{tf) — xj. 

For most path-integrals in physics, the path integral is 
quadratic in momenta and because of that it is customary to 
integrate first the momenta p{t) to obtain the Lagrangean of 
the problem, and be left with the path integral over the coor- 
dinate x{t). This route can be taken here too, but the resulting 
Lagrangean will have a position dependent mass, which un- 
necessary complicates the remaining integration over the co- 
ordinate. It is important to notice that the problem is linear 
in the coordinate x{t) and quadratic in p{t), so the original 
Schrodinger formulation will be simpler in the momentum 
representation. In the path-integral formulation, it turns out 
that the integration over x{t) is trivial and gives a delta func- 
tional. 



Both the continuous expression (|94} and the discrete one 
(I92> have their own computational advantages. In the follow- 
ing, I will use the discrete form ( I92L because for this problem, 
it allows for a quick and transparent solution. The exponent 
in \92\ can be rearranged in the form 



M 



S = i(pMXf -pQXi) - iejd^pj 



(96) 



M-l 



{Pj-i -Pj) +e-lPj + 



where I made the notation po — Pi ~ e(7Pi + iK^Pi/i). 

The path integral over the coordinate can be taken one by 
one, starting with the one over dxM-i- The result is a delta 
function 2n5{pM-i — Pm + (-IPm + ^c'^'^Pm/^)- The inte- 
gration over pm-1 resolves the delta function for the value 
of pM_i equal to pM-i = (1 — eT)pM, where the operator 
T = 7(-) + fc^(-)^/2. The integral over dxM-2 is done next, 
with the result 27r^(pM-2 — (1— followed by the inte- 
gral over dpM-2- The procedure is repeated until the last pair 
dxidpi. In this way we get po — {^ — £T)pi = (1 — eT)^^pM- 

The value of the transition probability is 



P{xf,tf \Xi,U) 



dpM i(pMXj-poXi)-ie'fe'^'^''^_^pj 



2tt 



(97) 

where in (|97}, pj = (1-ef for j = 0, M. Because 

of the approximation used in ( I90> , Eq.(|^ is valid only when 
the number of time splittings M oo. In this limit we have 



Pj - Pj-i 



Tpj 



dpt 
dt 



fpt with p{tf) 



The differential equation dpt /dt — Tpt 
Bernoulli equation and has the solution 

1 



Pt 



Vujr ^ 27 2'V 



= Pm 
(98) 

IPt + i^Pt is a 



(99) 
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The expression for the transition probability in the 
M oo limit is 



P{Xf,tf I Xi,ti) = 



dpM i(PMX I -poXi)-i^0 J*^ dtpt 



27r 



(100) 

where pt is given by Eq.j99ll. The time integral in JlOOt can 
be evaluated as 



dtpt 



2i 

1^1 



dy 



(ey - 1 



-7T 



C-1 



where 



2ij 



and T = t 



/ 



(101) 



(102) 



After minor rearrangements in the argument of the logarithm, 
the final result is given in the form of a Fourier integral 



P{xf,tf I Xi,ti) 



dp 



M 



2tt 



i{PM 



Xf-pQXi 
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FIG. 30: The complex plane of p. The branch cut due to the denomi- 
nator of Eq. <106t is shown as a broken line. The contour C has been 
deformed so that it encircles the branch cut. 



X exp 



276* 



In 



27 



-7T 



)P 



M 



where 



Po = p{t = U) 



PAie 



(103) 



(104) 



M 



Although, the integral from Ea.( ll03> seems complicated, it 
can be evaluated exactly by integration in the complex plane. 
Introduce the following notations 



The integral in ( I103> becomes 

r+oo cxp \ivxf 



and A = 27/k2(1 - exp(-7r)). (105) 

} 



27r 



(1 + ivi\y 



(106) 

There is a branch point singularity in the complex plane of 
p which is situated on the imaginary axis, for -p — iX. We 
take the branch cut to extend from iX to ioo. For positive 
Xf, the contour of integration is deformed so that it encircles 
the branch cut. The complex plane of p, and the contour of 
integration are presented in FIG|30| The contribution from 
the contour C can be split into two parts, one from the part 
to the left of the cut C_ and the other part, to the right of the 
cut C+. The Bessel function J^{z) has the following integral 
representation due to Schlafli 

where the integral is taken along a contour around the branch 
cut on the negative real axis, encircling it in the counterclock- 
wise direction. By making the change of variable p —ip in 
( I106> . the transition probability 

P{xf,tf \xi,ti) can be expressed in terms of the Bessel func- 
tion as 



P{xf,tf \xi,ti) X^e' 



(108) 



^i{2i^x^XfX'^ exp(-7T)). 



The modified Bessel function I,^{z) is defined in terms of the 
usual Bessel function Ji,(z) of imaginary argument by the re- 
lation 

I^{z) = e-T^ Jj.(e"/2^) for - tt < arg{z) < tt/2. 

(109) 

This relation allows to write the final result as 



P{xf,tf\x,,U) = X''e~^^^f+^''"'"^ 



(110) 



X I — A^e-^^j /,_i(2y'a;,a;/A2exp(~7T)), 

were A and ly are given by ( I105> . and T ~ tj — ti. 
The normalization of the probability density. 



dxf P{xf,tf I Xi,U) = 1, 



can be proved easily by making use of the result 



^ dxx-^+'e-'^^ Lm = j^ 



(111) 



We can check the validity of expression ( II lOl l by taking 
several limiting cases: 

• In the limit T 00, the transition probability 
P{xf,tf \xi,ti) should become the stationary probability dis- 
tribution (given in our paper). The argument of the Bessel 
function from il 10> goes to zero in the limit T — > 00. The 
modified Bessel function I^{z) has the following small argu- 
ment behavior 



for z 0^ 



r(M + i) 

Substituting the expression (fTT2l into ( fnol . we obtain 



(112) 



A'' 



P{xf,tf I a;,, -00) = P^ixf) = f^T-^Xf 



"^h-^'^f. (113) 



which is the stationary distribution. 

• In the short time limit, T (A — > 00), the 
probability distribution should approach a delta function, 
P{xf,ti I Xi, ti) — 6{xf — Xi). The argument of the Bessel 
function from ( II 10> goes to infinity in the limit T — > 0. The 
modified Bessel function (z) has the following large argu- 
ment behavior 



/27rz 



for z ^ 00 



(114) 



Substituting dl 141 back into dl lOt . we obtain 



P{xf,tf\xi,U) = f{xi/xf) lim 



A^oo V 47ra: 



S{xf - Xi) 



(115) 
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